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I. 


INTRODUCTION 


Finite  element  models  are  used  to  predict  the  dynamic  response  of  the  system  to 
various  excitations,  boundary  conditions  and  parameter  changes.  When  tests  are 
performed  to  validate  the  FE  model,  inevitably  their  results,  commonly  natural 
frequencies  and  mode  shapes,  do  not  coincide  with  the  expected  results  from  the  FE 
model.  These  discrepancies  stem  from  uncertainties  regarding  the  governing  equation  of 
the  system,  simplifying  assumptions  in  the  derivation,  and  inaccurate  boundary 
conditions.  Clearly  one  would  like  to  have  a  better  model,  based  on  both  the  theoretical 
and  the  experimental  results.  The  problem  of  how  to  modify  the  theoretic  model  from  the 
experimental  results  is  known  as  model  updating  (Kenigsbuch  and  Halevi,  1998). 

Finite  Element  Models  are  defined  by  a  large  number  of  physical  parameters,  but 
a  modal  test  of  a  structure  results  in  a  small  number  of  modal  parameters  which  are  used 
to  modify  the  model  parameters.  That  is  why  updating  a  finite  element  model  in  order  to 
produce  a  reliable,  confident  prediction  of  the  structural  response  of  a  system  tends  to 
become  a  difficult  task  to  be  accomplished. 

A  traditional  problem  associated  with  updating  a  FEM  has  been  obtaining 
sufficient  measurements  in  order  to  obtain  accurate  enough  results.  In  theory  it  should  be 
possible  to  measure  as  many  mode  shapes  as  required.  This  is  not  feasible,  because  the 
available  transducers  and  data  acquisition  hardware  limit  the  frequency  range  that  can  be 
measured.  Also,  as  the  frequency  increases  so  does  the  modal  density,  causing 
considerable  difficulty  for  the  modal  extraction  algorithms.  Even  if  these  difficulties 
could  be  overcome  it  would  still  be  impossible  to  accurately  measure  the  infinite  number 
of  modes  of  a  structure.  A  finite  element  model  only  accurately  predicts  the  lower  third  or 
so  of  the  natural  frequencies  and  mode  shapes.  Therefore  the  structure  could  not  be 
expected  to  reproduce  all  the  modes  predicted  by  a  model  (Ewins,  2000). 

FEM  can  be  very  large,  extending  in  many  cases  to  several  hundred  thousand 
Degrees  of  Freedom  (DOF),  or  more.  Experimental  modal  analysis  rarely  uses  more  than 
a  couple  of  hundred  transducers;  consequently  not  all  the  degrees  of  freedom  in  the 
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analytical  model  will  be  measured.  There  are  often  internal  nodes  which  cannot  be 
measured.  There  are  two  ways  of  overcoming  this  difficulty,  namely  reducing  the 
analytical  model  order  or  expanding  the  measured  mode  shapes. 

Procedures  have  been  developed  for  measuring  larger  experimental  databases. 
One  method  is  to  identify  additional  and  distinct  mode  frequencies  from  the  same  modal 
test  that  was  already  performed  to  identify  the  mode  frequencies  of  the  standard  system, 
without  the  need  for  physical  modification  of  the  structure.  The  additional  frequencies 
and  mode  shapes  correspond  exactly  to  the  mode  frequencies  found  when  combinations 
of  measured  coordinates  are  constrained  to  ground.  These  additional  frequencies  are 
therefore  associated  with  different  boundary  conditions  for  the  structure  and  no  physical 
change  in  the  boundary  conditions  has  been  made.  (Gordis,  1999). 

Since  the  boundary  conditions  are  not  actually  physically  applied,  they  are  termed 
Artificial  Boundary  Conditions  (ABC).  The  ABC  are  imposed  on  the  FEM  and 
correspond  to  ideal  constraints.  With  only  one  experimental  database,  this  design  yields  a 
separate  FEM  for  each  ABC  configuration.  With  the  exception  of  the  imposed  boundary 
conditions,  each  model  is  identical. 

In  a  set  of  spatially  incomplete  Frequency  Response  Function  (FRF)  data,  the 
ABC  are  the  constraints  that  define  an  Omitted  Coordinate  System  (O-SET).  In 
performing  a  vibration  test,  a  choice  is  made  as  to  the  set  of  coordinates  to  instrument 
with  response  transducers.  This  set  of  DOF  correspond  to  the  coordinates  that  are  going 
to  be  actually  measured. 

Sensitivity  based  updating  makes  use  of  the  modal  parameters  to  help  locate  a 
difference  between  a  FEM  and  an  existing  structure.  The  parameters  experimentally 
determined  from  the  actual  structure  are  the  natural  frequencies  of  the  structure,  the 
corresponding  mode  shapes,  and  the  modal  damping  ratios;  the  parameters  needed  from 
the  initial  FE  model  are  the  natural  frequencies  and  the  associated  frequency  sensitivities. 
Finite  element  model  sensitivities  are  the  first  derivative  of  the  eigen-problem  solution  set 
with  respect  to  the  design  variables  under  consideration. 
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Frequency  sensitivities  are  the  basis  for  a  linear  approximation  to  compute  the 
change  in  the  natural  frequencies  of  a  model  based  on  a  change  in  a  given  design 
variable.  Trying  to  locate  the  source  of  this  difference,  i.e.  either  the  error  in  the  model  or 
damage  in  the  structure,  is  a  problem  that  is  generally  underdetermined  due  to  the  number 
of  parameters  that  can  be  altered  and  the  limited  number  of  measured  parameters.  With 
the  application  of  ABC’s  a  better  solution  to  this  problem  can  often  be  formulated.  This 
method  of  frequency  updating  greatly  reduces  the  computational  time  required  to 
calculate  the  results  of  a  change  to  large  finite  element  models  by  eliminating  the  need  to 
resolve  the  eigen-problem  for  each  iteration.  This  process  has  been  used  extensively  for 
model  updating  (see  for  example,  Flanigan,  1998,  Luber,  1995,  Dascotte,  1995). 

Furthermore,  there  are  some  other  approaches  to  search  within  the  solution 
domain  for  potential  solutions,  such  as  the  Constrained/  Unconstrained  Optimization 
Routines.  The  Optimization  Routine  employed  searches  for  the  minimum  value  of  an 
objective  function  which  is  constructed  to  minimize  the  differences  between  the 
analytical  model  parameters  and  the  experimental  data.  This  organized  approach  of 
searching  for  the  optimal  combination  of  design  variables,  minimizes  the  number  of 
combinations  that  are  investigated  thereby  reducing  the  computational  requirements  of 
this  process. 

The  real  difficulty  to  this  method  of  correcting  a  finite  element  model  is  to 
determine  which  candidate  solution  most  closely  matches  the  experimental  data.  Multiple 
combinations  of  design  variable  changes  can  result  in  the  same  changes  to  the  natural 
frequencies  of  the  model  so  the  accuracy  of  the  frequencies  is  not  the  only  evaluation 
factor. 

An  improved  method  of  solution  evaluation  is  to  determine  the  effects  on  the 
analytical  mode  shapes  and  compare  them  to  the  experimental  mode  shapes.  This 
comparison  of  mode  shapes  is  known  as  the  Modal  Assurance  Criterion  (MAC)  and  is 
discussed  in  (Allemang,  1982).  The  MAC  is  a  measure  of  how  closely  the  analytical  and 
the  experimental  mode  shapes  coincide.  The  MAC  has  a  value  from  0  to  1,  where  the 
candidate  solution  with  a  MAC  value  closer  to  one  when  compared  to  the  experimental 

data  most  closely  resembles  the  experimental  case  and  is  selected  as  the  optimal  solution. 
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The  more  meticulous  the  Objective  Function  is  assembled,  that  is,  the  more 
properties  of  the  structure  are  considered  (design  variables),  the  better  the  error 
predictions  are  to  be  found  in  the  Constrained  Optimization  Routine. 

Three  different  Objective  Functions  were  evaluated  in  the  present  thesis  in  order 
to  correlate  the  confidence  in  the  damage  detection  results. 
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II.  THEORY 


The  matrix  equation  of  motion  for  a  spring  mass  oscillator  (SMO)  will  represent 
the  starting  point  of  the  derivation: 

[M]{x}  +  [C]{x]  +  [£]{*}  =  {F(t)}  (2.1) 

Excluding  the  damping[C] ,  equation  (2.1)  becomes: 

[M]{x}  +  [^]{x}  =  {F(0}  (2.2) 

Rearranging  equation  (2.2)  in  a  more  explicit  form,  it  is  written  as: 
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where  [M]  and  [K]  represent  the  mass  and  stiffness  matrices  which  define  the  physical 
properties  of  the  structure.  Both  matrices  are  square  and  symmetric,  that  is  of  size  n  x  n 
where  n  is  the  number  of  DOF  in  the  model.  The  vector  |F(t)}  represents  the  excitation 
vector  of  size  nxl  along  the  x  coordinate. 

Any  continuous  system  has  a  continuous  arbitrary  distribution  of  mass  and 
stiffness;  hence  it  also  comprises  an  infinite  number  of  DOF.  During  modal  tests,  only  a 
finite  number  of  DOF  are  measured,  and  these  measured  coordinates  represent  the 
analytical  set  or  (ASET),  where  the  measuring  instruments  (accelerometers)  are 
positioned.  The  remaining  unmeasured  DOF  are  considered  the  omitted  coordinate  set  or 
(OSET). 
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A. 


SPATIALLY  COMPLETE  AND  INCOMPLETE  DATA 


It  is  not  physically  possible  to  develop  a  complete  finite  element  model  of  a 
structural  system  over  a  frequency  domain.  This  would  require  a  number  of  measuring 
devices  (accelerometers)  to  be  equal  to  the  number  of  DOF  of  the  finite  element  model, 
hence  generating  it  to  have  a  determined  problem  where  the  numbers  of  equations  is 
equal  to  the  number  of  unknowns,  therefore  dealing  with  a  set  of  spatially  complete  data 
to  solve  the  problem  with  reliable  accuracy.  However,  practical  realities  only  allow  a 
limited  number  of  DOF  to  be  measured  in  a  modal  test.  Hence  the  result  for  a  real  world 
structure  is  a  measured  Frequency  Response  Function  (FRF)  matrix  that  is  spatially 
incomplete,  involving  fewer  measured  DOF  than  model  DOF,  therefore  obtaining  an 
underdetermined  problem. 

For  this  reason,  to  perform  the  structural  system  analysis,  that  is  either  the  model 
updating  or  the  damage  detection,  the  analytical  mode  shapes  of  the  Frequency  Response 
Function  matrix  must  be  reduced  in  size  such  that  only  the  mode  shapes  identified  from 
the  test/  measured  Frequency  Response  Function  matrix  are  the  ones  employed  to 
perform  the  system  analysis;  this  is  done  via  a  Reduced  Order  Model. 

The  Reduced  Order  Model  makes  use  of  the  ASET  and  OSET  coordinates,  where 
as  stated  before,  the  OSET  coordinates  are  not  associated  with  the  measurement  locations 
on  the  structure,  and  the  reduced  Frequency  Response  Function  matrix  is  composed  of 
response  information  from  the  limited  number  of  measurements  corresponding  to  the 
ASET  coordinates.  It  is  important  to  note  that  the  reduced  FRF  matrix  implicitly  defines 
a  dynamically  reduced  Impedance  matrix,  furthermore,  because  the  basic  system  has  not 
been  altered  and  that  it  still  has  the  same  number  of  DOF,  even  though  it  was  determined 
not  to  describe  all  of  the  DOF,  the  elements  remaining  in  the  reduced  FRF  matrix  are 
identical  to  the  corresponding  elements  of  the  full  FRF  matrix.  (Ewins,  1984).  The  matrix 
still  contains  all  the  modal  information  of  a  fully  described  matrix. 
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B. 


ANALYTICAL  AND  OMITTED  COORDINATE  SETS 


Rearranging  equation  (2.3)  into  the  frequency  domain  for  the  steady-state 
harmonic  response  for  a  full  order  model  without  the  damping,  it  becomes: 
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where  Q  is  the  frequency  of  the  harmonic  excitation.  Considering  that  the  number  of 
DOF’s  are  comprised  for  the  measured  and  unmeasured  set  of  coordinates,  that  is  the 
ASET  and  OSET  coordinates,  we  can  express  equation  (2.4)  in  the  following  form: 
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(2.5) 


where  the  subscript  “a”  refers  to  the  ASET  coordinate  and  “o”  to  the  OSET  coordinate. 
Equation  (2.5)  is  also  known  as  the  impedance  model  and  it  can  be  redefined  as: 


where  [Z]  represents  the  system  impedance  matrix  evaluated  at  the  forcing  frequency,  Q . 
Assuming  a  static  constraint,  that  is  no  excitation  on  the  omitted  coordinates 
then{/0}  =  {0} ,  and  equation  (2.5)  is  rewritten  as: 
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Rearranging  equation  (2.7)  into  two  separate  equations  and  substituting  the  OSET 
coordinates  in  terms  of  the  ASET  coordinates,  only  one  equation  is  required: 

K,  ]R }+ K,  T  „  1-  sv  ]k }+ [m„  ]{*„ }] = o  (2.8) 
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Solving  equation  (2.8)  in  terms  of  {x0}  it  becomes: 


{x  }  =  \i  -  n2K-'M00 1"1  [-  K~lKna  +  Q 2K-lMoa  \xa 
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(2.9) 


or 


W=KT'K]W 


(2.10) 


By  definition  the  inverse  of  the  OSET  impedance  matrix  [ZGO]  can  be  expressed  as: 


[/  -  Q 2K-'Mnn  I  '  = - r - l- - - - q  Adj\l  -  Q.2K;\Mnn  1  (2.11) 

L  oo]  Det[I  -Q.2 K;'M JL  oo  oo] 

where  the  Det[*  ]  indicates  the  determinant  and  Adj[*  ]  indicates  the  adjoint  matrix. 
From  equation  (2.11),  it  is  noticeable  that  the  bracketed  inverse  term  does  not  exist  at 
those  frequencies  which  satisfy: 

D«[/-n2Oc]=o  (2.12) 


By  definition  the  eigenvalues  which  satisfy  equation  (2.12)  are  those  defined  by 
the  [A^]  and  \M00\  matrices,  hence  the  resulting  frequencies  correspond  to  the  OSET 
frequencies  only. 

Considering  that  the  eigenvalues  and  eigenvectors  of  a  system  are  defined  by  the 
unrestrained  DOF  of  the  corresponding  system,  the  OSET  coordinates  of  the  OSET 
system  are  unrestrained  and  the  ASET  coordinates  are  constrained  to  ground. 

The  FRF  matrix  derivation,  as  well  as  the  Reduced  Order  Model  are  described  in 
the  following  sections. 


c. 


REDUCED  ORDER  MODEL 


The  process  of  instrumenting  a  structure  with  a  finite  number  of  response 
transducers  defines  a  Reduced  Order  Model,  where  the  impedance  of  the  Reduced  Order 
Model  is  nonlinearly  dependent  on  the  impedance  of  the  Full  Order  Model  (Gordis, 
1996). 

In  order  to  perform  model  updating  using  spatially  incomplete  data,  a  reduction  of 
the  model  data  is  done,  that  is  the  reduction  of  the  analytical  model,  with  no  loss  of 
generality  of  the  results. 

Considering  a  full  order  “exact”  FRF  model  of  a  structure,  we  would  have  a  FRF 
matrix  of  infinite  dimensions: 
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where  the  experimental  FRF  matrix  actually  measured  in  a  real  test  is  seen  to  be  a  matrix 
partition  which  has  been  extracted  from  the  infinite  dimension  matrix,  i.e. 


[tfj 


(2.14) 


The  overbar  in  equation  (2.14)  specifies  a  reduced  model,  the  superscript  x 
denotes  an  experimentally  measured  FRF,  and  the  \Haa]  matrix  represents  a  structural 

dynamic  model  that  has  been  reduced  using  exact  dynamic  reduction  (Gordis,  1996). 
From  the  partitioned  inverse  relation  of  the  FRF  matrix  to  its  associate  impedance, 
i.e.[Zl[tfl  =  / ,  we  have: 
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Now,  solving  for  [Haa]  gives: 


[hJ  =[z. 


(2.15) 
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(2.16) 


Since  it  was  assumed  no  excitation  on  the  omitted  coordinates  ({/o}  =  {0}),  the 

exact  dynamic  relationship  between  the  OSET  coordinates  and  the  ASET  coordinates 
gives: 


Rearranging  equation  (2.5)  using  equation  (2.17),  the  dynamic  reduced  model 

results: 

if)=\z  -Z  Z_1Z  1{jc  }  (2.18) 

Comparing  equations  (2.16)  and  (2.18),  we  can  see  that  the  partition  of  the  infinite 
dimensional  FRF  matrix  measured  in  a  test  is  identically  equivalent  to  the  inverse  of  a 
dynamically  reduced  impedance  matrix.  Again  we  see  the  presence  of  the  OSET  system 
dynamics  in  the  termZ^1.  The  formula  of  the  matrix  inverse,  given  by  equation  (2.11) 

applies  in  the  same  fashion,  and  since  every  element  in  the  Zj  is  singular  at  the  natural 
frequencies  of  the  OSET,  it  is  seen  from  equation  (2.16)  that  the  elements  of  H~fa  will  be 
singular  at  the  OSET  natural  frequencies  (Gordis  1999). 

Given  that  the  measured  FRF  matrix  implicitly  defines  a  dynamically  reduced 
impedance  model,  we  pursue  the  reduction  of  the  FEM  in  the  same  manner.  That  is,  we 
calculate  an  nxn  FRF  matrix  from  the  FEM,  and  simply  extract  the  portion  of  this  matrix 
which  corresponds  to  the  coordinates  measured  in  a  test.  Hence  the  resulting  extraction- 
reduced  finite  element  FRF  matrix  retains  all  the  modal  content  of  the  original  model. 

D.  FREQUENCY  RESPONSE  FUNCTION  [H] 

The  Frequency  Response  Function  is  a  characteristic  of  the  system  where  the 
frequency  response  is  the  frequency  domain  input-output  relationship  which  can  be 
considered  the  equivalent  of  the  time  domain  impulse  response  function  (Ewins,  2000). 
The  FRF  function  contains  the  magnitude  and  phase  (complex  amplitude)  of  the  response 

at  the  corresponding  DOF  “z  ”  due  to  a  harmonic  force  of  unit  amplitude  at  DOF  “j  The 
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time  domain  measured  data  is  transformed  to  the  frequency  domain  using  the  Fast 
Fourier  Transform  (FFT)  algorithm. 


The  inverse  of  the  FRF  matrix  ur  is  known  as  the  impedance  matrix  [Z] ,  that 


is: 


[Zip)}  =  [HP)]-1 

and  as  it  was  implied  from  equations  (2.5)  and  (2.6): 

[Z(  Q)]  =[K-Cl2M  +  j€lC\ 
Analyzing  equation  (2.6)  in  terms  of  an  n  DOF  system 
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Transferring  from  physical  to  modal  coordinates  let  {x}  =  [d>]{g} ,  where  O  is 
the  mass  normalized  mode  shapes  of  the  system  and  g  represents  the  generalized 
coordinate  set,  and  rewriting  equation  (2.21) 

[Z][d)fo}={/}  (2.22) 

Pre-multiply  by  [O ] 7  and  expanding  [Z]  from  equation  (2.20) 

[[®f  [X][®]-fi2  [®f  [M][4>]+;n[4>f  [C][®]]{?)  =[®f  {F}  (2.23) 

Using  orthogonality  [o]r  [M  ][o]  =  1 ,  and  assuming  proportional  damping 
[C]  =  cc*[K]  +  j3*[M]: 

[®,2  -n2  +2j$l<o^{q}  =[®f  {F}  (2.24) 
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Where  &>.  is  the  natural  frequency  of  the  ith  mode,  C,  is  the  damping  ratio  and 
[a>f  -Q2  +2j^D.(oi~^  is  known  as  the  modal  impedance  matrix  and  is  diagonal.  This 
matrix  is  inverted  to  find  the  modal  frequency  response: 


M=M 


[&>; -n’ +2/rriffil] 


(2.25) 


Transforming  back  equation  (2.25)  into  physical  coordinates  by  pre-multiplying 
by  [<i>]  and  using  {3}  =  [0]r  {F} ,  the  modal  decomposition  of  the  FRF  matrix  in 
physical  coordinates  is  as  follows: 


M  =  [®] 


[cof  -Q.2  +2jCno)l~\ 
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(2.26) 


where  [H (Q)]  can  be  defined  as: 


[#(«)]  =  Nr  2  n2  T  —  1 

y(oj  -Q  +  2j£Q.coi  \ 

Writing  equation  (2.27)  in  summation  form: 
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Or  for  any  element  in  terms  of  [//(>.(Q)J  : 


(2.28) 
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E.  DRIVING  POINT  AND  TRANSFER  POINT  FUNCTIONS 

The  Frequency  Response  Function  represents  a  harmonic  response  in  the 
ith  coordinate,  caused  by  a  single  harmonic  force  applied  either  at  the  same  ith  coordinate 
or  at  a  different  j,h  coordinate.  Hence  the  Driving  point  Function  is  one  where  the 
response  coordinate  and  the  excitation  coordinate  are  identical;  where  as  the  Transfer 
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Function  is  one  where  the  response  and  excitation  coordinates  are  located  at  different 
DOF  points.  The  FRF  is  the  foundation  of  the  modal  testing.  It  shows  a  direct  connection 
between  the  modal  properties  of  a  system  and  its  response  characteristics;  hence  it 
provides  an  efficient  means  of  predicting  responses  (Ewins,  1984). 

The  following  example  is  a  two  DOF  system  with  no  damping  (Ewins,  2000).  It 
displays  the  characteristics  of  both  a  driving  point  and  transfer  function  for  the  complete 
FRF. 


Assembling  the  stiffness  matrix: 

Kx+K2  -k2 

-k2  k2+k3 

The  mass  matrix  is  defined  as: 

AT,  0  " 

0  M2 


For  the  purpose  of  this  example  MX=M2  =  1.0  and  Kx  =  K3  =  0.4  and  if2=0.8. 
This  yields  the  following  results: 


Natural  frequencies:  {< co(rad  /  sec)} 
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mod  es 
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k= i  o)k  — q  +  2jCfiojk 


For  the  Driving  Point  function 
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For  the  Transfer  function. 


[Ff12(Q)] 
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It  is  interesting  to  determine  what  controls  whether  a  particular  FRF  will  have 
positive  or  negative  modal  constants  and  thus  whether  it  will  exhibit  antiresonances  or 
not.  As  it  was  noticed  from  equation  (2.29),  it  is  the  product  of  two  eigenvectors 
elements,  one  at  the  response  point  and  the  other  at  the  excitation  point.  Clearly,  if  we  are 
to  consider  a  driving  point,  then  the  modal  constant  for  every  mode  must  be  positive,  it 
being  the  square  of  a  number.  This  means  that  for  a  point  FRF,  there  must  be  an 
antiresonance  following  every  resonance,  without  exception.  This  is  displayed  in  the  top 
portion  of  Figure  2  which  is  the  Driving  Point  FRF[/7n] . 


The  situation  for  the  Transfer  Function  is  somewhat  different  because  clearly  the 
modal  constant  will  sometimes  be  positive  and  sometimes  negative,  depending  upon 
whether  the  excitation  and  response  move  in  phase  with  each  other  or  not.  Thus  we 
expect  Transfer  FRF  measurements  to  show  a  mixture  of  antiresonances  and  minima. 
However,  this  mixture  can  be  anticipated  to  some  extent  because  it  can  be  shown  that,  in 
general,  the  further  apart  are  the  two  points  in  question,  the  more  likely  are  the  two 
eigenvector  elements  to  alternate  in  sign  as  one  progresses  through  the  modes  (Ewins, 
2000).This  is  shown  in  the  lower  plot  of  Figure  2  which  is  the  H12.  Conversely  in  the 
region  between  the  two  natural  frequencies  the  two  terms  are  additive  and  thus  do  not 
create  an  antiresonance. 
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Driving  Point  FRF,  Fill 


Transfer  FRF,  HI 2 


Figure  2.  Two-DOF  Frequency  Response  Function. 

The  principles  demonstrated  in  the  above  example  can  be  applied  to  any  number 
of  DOF.  This  concept  of  the  existence  of  an  antiresonance  between  any  two  modes  of  a 
driving  point  FRF  is  of  great  significance  in  the  use  of  Artificial  Boundary  Conditions  as 
it  is  shown  in  the  next  chapter. 
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III.  ARTIFICIAL  BOUNDARY  CONDITIONS 
CONFIGURATION  OF  NATURAL  FREQUENCIES 


The  improvement  of  a  FEM  is  often  a  necessary  step  in  order  for  the  model  to  be 
used  with  confidence  in  the  prediction  of  structural  response.  This  improvement  is 
difficult  to  achieve  due  to  the  spatially  incompleteness  of  the  measured  FRF  matrix 
giving  an  underdetermined  problem.  Undoubtedly  there  exists  a  recognized  need  to 
increase  the  size  of  the  known  parameter  data  base;  it  has  been  found  that  a  number  of 
additional  mode  frequencies  can  be  identified  from  the  same  modal  test  without  the  need 
for  any  physical  modification  of  the  structure.  These  additional  and  distinct  mode 
frequencies  correspond  exactly  to  those  mode  frequencies  found  when  combinations  of 
measured  coordinates  are  restrained  to  ground.  This  last  idea  is  referred  as  the  ABC 
concept,  which  as  it  will  be  shown,  improves  the  updating  process  by  associating  the 
natural  frequencies  with  different  boundary  conditions. 


A.  RELATION  BETWEEN  THE  DRIVING  POINT  ANTIRESONANCES  TO 
ABC  FRQUENCIES 


Evoking  the  Driving  Point  Function  concept  discussed  in  the  last  chapter  it  can  be 
shown  that  the  ABC’s  frequencies  of  a  given  ASET  coordinate  are  defined  by  the 
corresponding  driving  point  antiresonances.  Furthermore,  the  ABC’s  mode  frequencies 
not  only  provide  a  greater  number  of  frequencies  for  the  system,  but  also  provide  a  means 
to  reduce  ill-conditioning  in  the  solution  of  the  sensitivity  equations. 

From  equation  (2.29)  the  driving  point  FRF  is 


modes 


k= 1 


(fr-)2 

col -Or 


(3.1) 


where  O,  is  the  mass  normalized  mode  shape  element,  oy-  is  the  kt h  natural  frequency, 
and  Q  is  the  forcing  frequency. 
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The  frequency  of  the  anti-resonance  of  HU(Q.)  is  given  by 


+  /?>,2 

*n+*n 


(3.2) 


where  the  modal  residue  is  giving  by  Ryk  =  {<!>(:,  k)}  { 0(: ,  k)}.  The  following  example 
(Gordis,  1999),  shows  how  frequency  is  identical  to  the  natural  frequency  of  the  system 
in  Figure  3,  to  the  driving  point  DOF  constrained  to  ground  as  shown  in  Figure  4. 


1.  Calculation  of  ABC  Frequencies  for  a  2-DOF  System 

Given  the  shown  undamped  system  of  Figure  3: 


Let  Mi  =  M2  =  1.0  and  Ki  =  K2  =  1.0,  the  frequency  of  the  single  antiresonance  of 
the  system  isQan<;>ej  =  V2  rad/sec,  which  is  identically  equal  to  the  single  natural 

frequency  of  the  ABC  system  of  Figure  4  which  co  =  \[l  rad/sec. 


Figure  4.  ASET,  DOF  #1  constrained  to  ground. 
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Mode  shapes:  {0*}  = 


Natural  frequencies:  {^(raJ/sec)} 
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F igure  5 .  Driving  Point  FRF ,  H 1 1 . 

The  Hn( Q)  plot  in  Figure  5  graphically  shows  the  calculated  mode  frequencies, 
as  well  as  the  anti-resonance  frequency  as  it  was  computed  from  equation  (3.2)  is 
Q antires  =  V2  rad/sec.  or  1.41  rad/sec. 

By  the  use  of  equation  (2.16),  the  frequencies  of  the  ABC  system  can  be  found 
just  by  calculating  inverse  of  //, ,  (Q) ,  where  the  ASET  is  pinned  at  DOF  #1,  and  H~'a  (Q) 


Driving  Point  FRF,  H11 
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is  the  inverse  of  HU(Q.) .  A  clear  correlation  is  depicted  in  the  bottom  graph  of  Figure  6, 

where  the  single  natural  frequency  corresponds  to  the  antiresonance  frequency. 
Furthermore,  this  singular  frequency  happens  to  be  the  natural  frequency  of  the  ABC 
system  obtained  by  constraining  DOF  #1  to  ground  (Figure  4).  Hence  it  is  shown  the 
power  and  usefulness  of  equation  (2.16),  where  the  driving  point  antiresonance 
correspond  directly  to  the  natural  frequencies  of  the  system,  with  the  driving  point 
coordinate  constrained  to  ground  (Gordis,  1999). 


•6 


o 


Plot  A.  Driving  Point  FRF,  H11 


rad/sec 


Plot  B.  Z  =inv(H1 1 )  of  the  ABC  System 


Figure  6.  Plot  A.  Driving  Point  HI  1(0)  of  system  1,  Plot  B.  [Haa(0)]‘'  of  system  2. 

2.  Calculation  of  ABC  Frequencies  for  a  Free-Free  Beam 

To  further  show  the  usefulness  of  the  ABC  concept,  a  free-free  beam  is  analyzed 
(Gordis,  1999).  The  present  model  depicted  in  Figure  7  is  examined,  it  consists  of  10 
beam  elements  and  each  element  contains  two  nodes  with  two  DOF  each,  giving  a  total 
of  22  DOF  system.  The  translational  DOF  correspond  to  the  odd  numbers  where  as  the 
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rotational  DOF  correspond  to  the  even  numbers.  Response  transducers  are  set  up  at  the 
translational  DOF:  1,  5,  9,  13,  17,  and  21  as  follows: 


Figure  7.  Transducers  located  at  DOF  1,  5,  9,  13,  17,  and  21. 

The  set  of  coordinates  chosen  to  be  measured  with  response  accelerators  represent 
the  ASET,  defined  as  [1  5  9  13  17  21].  It  can  be  assumed  that  the  excitation  is  applied  at 
each  of  the  ASET  DOF,  resulting  in  a  6  x  6  FRF  matrix.  The  impedance  matrix  H~'a  (Q) , 
the  scalar  inverse  of  Hn{ Q),  is  calculated  at  excitation  frequencies  from  0-800  hertz. 
Figure  8  shows  the  driving  point  FRF  of  the  system. 


Figure  8.  Driving  Point  FRF  ASET  DOF:  1  5  9  13  17  21.  [From  Gordis,  1999]. 
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The  ABC  system  frequencies  are  found  by  calculating  the  impedance  matrix 
H~al{ Q).  By  plotting  the  1,1  element  //,  ,(Q)  of  the  matrix,  the  ABC  system  frequencies 
can  be  seen  in  Figure  9: 


The  ABC  frequencies  correspond  exactly  to  the  natural  frequencies  of  the  system 
with  all  of  the  measured  coordinates  ASET  constrained  to  ground.  This  configuration  is 
shown  in  Figure  10: 


Figure  10.  ABC  Configuration  ASET  [1  5  9  13  17  21]  (Measured  coordinates 

restrained  to  ground). 
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IV.  SENSITIVITY  BASED  FEM  UPDATING  USING  ABC 


The  disparity  in  the  number  of  known  parameters  (measured)  versus  the  number 
of  parameters  to  be  adjusted  in  order  to  update  a  FEM,  defines  a  spatially  incomplete 
structure  as  well  as  an  underdetermined  problem. 

The  importance  of  obtaining  sensitivities  of  the  dynamic  characteristics  of  a 
system  with  respect  to  the  changes  of  its  parameters  lies  in  the  fact  that  they  are  critical 
for  achieving  efficient  design  modification,  and  for  predicting  the  optimal  structural 
modifications  within  a  prescribed  dynamic  characteristics  change.  In  such  a  case,  there 
will  be  numerous  possible  modifications  that  can  accomplish  the  change.  Sensitivity 
analysis  can  provide  the  most  effective  answer. 

From  a  mathematical  point  of  view,  sensitivity  analysis  implies  differentiation.  It 
can  be  carried  out  on  the  modal  data  from  the  FRF  of  a  dynamic  system.  The  variables  in 
this  analysis  are  usually  system  parameters  such  as  mass  or  stiffness  parameters.  The 
functions  of  the  analysis  can  be  natural  frequencies  and  mode  shapes  of  the  system 
(modal  data),  or  its  FRF.  The  assumption  of  infinitesimal  changes  demanded  by 
differentiation  applies  to  sensitivity  analysis.  Hence  the  sensitivity  based  updating 
concept  is  used  for  error  localization  in  order  to  find  the  parameters  requiring 
adjustments.  The  ABC  concept  can  be  used  in  addition  to  the  standard  baseline  system 
mode  frequencies  in  model  updating  and  damage  detection.  As  stated  previously,  the 
ABC  frequencies  correspond  to  the  same  structural  system  as  the  baseline  frequencies, 
but  with  different  boundary  conditions.  The  governing  equation  for  sensitivity  based 
updating  is: 

{a<u2}  =[T]{ADv}  (4.1) 

where  |Ao2|  is  a  vector  of  natural  frequency  errors.  The  errors  are  the  difference 
between  the  experimental  and  the  analytical  natural  frequencies  {<»*  -coa} .  The  {ADv} 
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term  is  the  vector  of  changes  to  be  calculated  for  the  specified  model  parameters,  known 
as  the  design  variables,  and  [ T ]  is  the  sensitivity  matrix. 

Each  column  of  the  sensitivity  matrix  represents  an  element  of  the  model,  where 
as  each  row  represents  a  mode  of  the  model. 

Each  ABC  system  defines  additional  rows  to  equation  (4.1)  that  is  it  adds  more 
modes,  defining  the  equation: 

|Aui2(,)|  =[r(,)]{ADv}  (4.2) 


where 


T(k) 

,J  dDVj 


and  cof(k)  represents  the  zth  natural  frequency  of  the  Mi  ABC 


configuration  system.  The  baseline  system  quantities  are  represented  with  the  superscript 
“0”  and  the  ABC  systems  are  identified  with  the  superscripts  from  1  to  “k”,  where  k 
represents  the  number  of  ABC  used  in  the  system.  Combining  baseline  and  ABC  system 
equations: 


Aco] 


2(0) 

’BASE 


A  co 


2(1) 


ABC  1 


A  co 


2  (k) 


ABCk 


f  t(°)  ^ 

1  BASE 


-(1) 


ABC  1 


T'ik) 

V  ABCk  J 


dV, 


dV. 


(4.3) 


The  degree  of  coupling  between  the  ABC  systems  and  the  baseline  system  in 
equation  (4.3)  can  be  adjusted  by  deleting  or  retaining  individual  columns  of  the  \j’k  J . 

Partial  coupling  between  the  baseline  system  and  the  ABC  system  is  established  if  the 
design  variables  (DV)  are  partitioned.  That  is,  some  of  the  DV  are  associated  only  with 
the  baseline  system,  {A DV0},  some  are  associated  with  both  the  baseline  and  the  ABC 
system,  {ADE0,1} ,  and  some  are  associated  only  with  the  ABC  system,  JADE1 } ,  (Gordis, 
1999).The  use  of  the  ABC  system  sensitivities  helps  to  eliminate  or  at  least  to  reduce  the 
problem  of  poorly  conditioned  or  rank-deficient  [T]  matrices.  Columns  of  [f0]  can  be 

replaced  with  columns  of  J  in  order  to  improve  the  conditioning.  Two  closely  spaced 
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elements  in  a  model  will  have  nearly  dependent  columns  of  [r0] ,  preventing  the 
discrimination  of  error  or  damage  between  two  elements.  One  column  of  [r0]  associated 

with  one  of  the  two  elements  can  be  replaced  with  a  column  of  [  Tk  J ,  and  the  associated 

baseline  system  natural  frequency  replaced  with  the  ABC  system  frequency  (Gordis, 
1999). 

A.  SENSITIVITY  MATRIX  DERIVATION  [T] 


Several  approaches  are  available  for  performing  a  sensitivity  analysis.  The 
following  derivation  is  based  on  the  parameterization  of  the  eigenvalue  problem  (Inman, 
1989).  Consider  an  n  DOF  system  defined  by: 

M(D  V)x(t)  +  K(D  V)x(t)  =  0  (4.4) 

where  the  corresponding  eigenvalues  problem  is: 


[A:-AiAf]{®/}  =  {0} 


(4.5) 


Differentiating  equation  (4.5)  with  respect  to  the  design  variable  {DV},  as  it  is  the 
vector  of  design  parameters  representing  the  change  in  mass  matrix  [M]  and/or  change 

in  stiffness  matrix  [ K] ,  which  are  adjustable  for  the  FE  model. 


A K 
A DV 


-A  i 


AM 
A DV 


Ah  ' 
A DV 


A  O,  1 


W  +  lK-Mt]  —  ={0} 


(4.6) 


Expanding  equation  (4.6)  and  pre-multiplying  by  {O  j  r  leads  to: 


{ch}7 


AK 
A DV 


{0,}-A{(Ej7 


AM 
A DV 


{®}- 


AA 
A DV 


And  invoking  the  matrix  identity  property  \a\T  [6]  {c{  =  |c}r  [6]{c/{  the  last  term 
on  the  left  hand  side  can  be  replaced  by: 


[  AO,  l7 

Wj 


[A  -/l,M]{Oi}  =  0 


(4.8) 
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Since  [ K  -  AM  ]  {0; }  =  {0}  the  overall  equation  is  reduced 


A K 
A DV 


AM 
A DV 


w- 


Ah 
A DV 


{Oi}3"  [M  ]  {O/}  =  {0}  (4.9) 


Using  orthogonality,  where  [<D;]r  [M][d>/]  =  1  equation  (4.9)  further  reduces  to 


A K 
A DV 


{<!>,}- A,{<D,}7 


AM 
A DV 


{«>,•}- 


Ah 

ADV 


=  M 


(4.10) 


Bringing  the  mass  and  stiffness  portion  to  the  right  hand  side 
Ah 


ADV 


-w 


AK  AM 
ADV  '  ADV 


{O,} 


(4.11) 


Finally  it  is  demonstrated  that  the  sensitivity  matrix  involves  differentiation  with 
respect  to  the  design  parameter  where  as  for  this  case  are  either  consider  mass  and/or 
stiffness,  and  the  function  of  the  analysis  is  with  respect  to  the  natural  frequencies  of  the 
system  as  it  is  shown  in  equations  (4.12)  and  (4.14): 

{<&/}  where  [A K]  =  [Kx-Ka]  (4.12) 

and 


Ah 

AK 

ADV 

=w 

ADV 

[-^s7(//neM  _  —  I®1'} 


AK 

ADV 


(4.13) 


and 


Ah 

ADV 


-h 


AM 

ADV 


where  [AM ]  =  [Me  -  Ma\ 


(4.14) 


\T  l  = 

L  mass  J  ( 


-h 


AM 

ADV 


{<!), 


(4.15) 
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B. 


BEAM  SPECIFICATIONS 


A  Finite  Element  computer  program  was  written  to  calculate  the  Sensitivity 
matrix  of  a  twenty  element  cantilever  beam  model,  which  consists  of  an  aluminum  bar, 
with  a  partially  drilled  and  tapped  10/32”  diameter  hole  at  the  free  end,  which  allows  the 
load  cell  to  be  attached  to  the  structure.  Table  1  provides  a  summary  of  the  beam 
specifications. 


PARAMETER 

VALUE 

Length 

42  inches 

Width 

1 .5  inches 

Thickness 

0.5  inches 

Density 

0.11  lb /in3 

Elasticity  Modulus 

10  E 6lbf  / sec2 -in  . 

Table  1.  Beam  Specifications. 


C.  MASS  AND  STIFFNESS  SENSITIVITY  ANALYSIS 

The  computer  program  built  to  assemble  the  sensitivity  matrix  calculates  the 
mode  frequency  sensitivities  according  to  equation  (4.11).  This  process  is  performed  in 
two  steps.  The  first  step  calculates  the  Mass  Sensitivity  matrix  and  the  second  calculates 
the  Stiffness  Sensitivity  matrix.  Hence  the  [T]  matrix  is  of  size  m  x  p,  m  being  the 
number  of  modes  and  p  being  the  number  of  elements. 

In  the  Mass  Sensitivity  calculation,  a  small  perturbation  of  1%  of  mass  was 
applied  to  each  element  of  the  beam;  this  is  done  using  equation  (4.14).  Each  subsequent 
column  of  [T]  was  calculated  at  every  element  of  the  model. 
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On  step  two  for  the  Stiffness  Sensitivity  matrix  calculation,  the  same  1%  stiffness 
perturbation  is  applied  iteratively  to  each  element  and  equation  (4.12)  is  used  to  assemble 
the  corresponding  part  of  the  Sensitivity  matrix. 

Once  the  two  steps  are  completed,  both  sensitivities  matrices  are  assembled 
[T(M)]  &  [T(K)],  hence  obtaining  a  [T]  total  of  size  m  x  2 p,  where  the  first  set  of  twenty 
columns  yield  the  changes  in  mass  and  the  second  set  o  twenty  columns  yields  the 
changes  in  stiffness. 

The  following  plots  show  a  normalized  mass  and  stiffness  sensitivity  distribution 
along  the  20  elements  of  the  beam  model.  Each  element  is  indicated  by  a  bar  which 
represents  its  corresponding  response  to  the  change  of  the  system  natural  frequency,  with 
respect  to  the  corresponding  change  of  the  design  variable. 

It  will  be  seen  in  an  overall  sense  on  all  of  the  Figures  below  how  significantly  the 
sensitivity  is  influenced  by  the  mode  shapes.  It  is  clear  from  Figure  11,  how  the  Base 
Stiffness  Sensitivity  is  more  likely  to  detect  a  change  in  the  natural  frequency  closer  to 
the  fixed  end  of  the  cantilever  beam  model,  rather  than  at  the  free  end,  where  the  red 
sensitivity  bars  are  almost  zero. 
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In  Figure  12,  it  is  evident  how  the  Base  Mass  Sensitivity  is  more  likely  to  detect  a 
change  in  the  natural  frequency,  closer  to  the  free  end  of  the  cantilever  beam  model, 
rather  than  at  the  fixed  end,  where  the  magenta  sensitivity  bars  are  almost  zero. 
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It  should  be  noted  that  the  Mass  Sensitivity  rows  are  represented  in  absolute 
values,  therefore  an  increase  of  mass  would  generally  lower  the  system  natural 
frequencies,  but,  it  was  preferred  to  show  the  influence  with  respect  to  element  position 
and  hence  all  magnitudes  were  made  positive. 

The  next  four  figures  show  a  particular  trend  followed  by  the  Stiffness  and  Mass 
Sensitivities  when  ABC’s  are  applied  to  the  system. 
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The  green  pointer  with  triangular  shape  at  element  twelve  of  Figures  13  and  15 
and  element  four  of  Figures  14  and  16,  represent  the  pinned  ABC  applied  at  DOF  25  and 
9,  respectively  as  follows: 


Figure  13.  Stiffness  Sensitivity  distribution  with  ABC  at  DOF  25. 


Figure  14.  Stiffness  Sensitivity  distribution  with  ABC  at  DOF  9. 
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Figure  15.  Mass  Sensitivity  distribution  with  ABC  at  DOF  25. 


Mode  1 


Mode  3 


Mode  5 


- 1 - 1 - 

a □ ^ □ 

- 1 - 

□ □ 

- 1 - 

^  n 

- 1 - 

n  .  1 

1 

1 - 

□ 

- 

□  2  4  6  3  10  12  14  16  13  20 


Figure  16.  Mass  Sensitivity  distribution  with  ABC  at  DOF  9. 
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In  general  when  ABC’s  are  applied  to  the  Stiffness  Sensitivity  matrix,  the  model 
has  a  tendency  to  increase  the  values  of  their  sensitivities  at  that  particular  “pinned” 
ABC,  whereas  for  the  Mass  Sensitivity,  the  pinned  ABC  seems  to  drive  the  sensitivity 
values  away  from  the  DOF’s  where  fixed  conditions  are  emulated;  these  fixed  conditions 
can  be  thought  as  those  representing  the  fixed  end  of  the  beam  as  well  as  the  pinned  ABC 
set  up  at  the  corresponding  DOF’s  25  and  9  as  it  was  showed  on  Figures  15  and  16, 
respectively. 

In  summary,  the  previous  Figures  expressed  that  in  places  of  low  sensitivity,  the 
beam  model  is  not  very  likely  to  predict  the  changes  in  natural  frequencies  due  to  the 
change  in  the  design  variable,  and  that  is,  in  a  real  world  scenario  those  low  sensitivity 
areas  would  be  the  places  where  damage,  if  any,  is  not  going  to  be  accurately  detected. 
Hence,  the  sensitivity  concept  leads  to  a  key  role  in  the  updating  of  a  FEM.  Furthermore, 
the  sensitivity  analysis  showed  us  two  important  points.  The  first  is  the  big  influence  of 
the  mode  shapes  into  the  sensitivity  distribution  of  the  model,  and  two;  how  the  correct 
manipulation  of  the  ABC’s  especially  with  regard  to  the  stiffness  sensitivity  portion, 
could  lead  to  a  better  prediction  of  the  errors  in  the  FEM  updating. 
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V.  KINETIC  AND  STRAIN  ENERGY  CORRELATION  TO 

SENSITIVITY  MATRIX 


Because  the  natural  frequencies  and  mode  shapes  of  a  structure  are  dependent  on 
the  mass  and  stiffness  distributions,  any  subsequent  changes  in  them  should, 
theoretically,  be  reflected  in  changes  in  the  frequency  and  mode  shapes  of  the  structure. 
Hence,  the  first  step  in  a  vibration  problem  is  to  determine  the  important  natural 
frequencies  of  the  system;  there  are  two  general  methods  to  determine  the  natural 
frequencies  of  the  system:  The  Newtonian  and  the  Energy  method  approach. 

The  Rayleigh’s  Energy  method  employs  an  energy  balance  which  is  essentially 
the  scope  of  this  chapter.  For  the  sake  of  the  following  development,  the  Euler-Bemoulli 
beam  theory  is  utilized,  where  the  Euler-Bemoulli  equation  for  beam  bending  is: 


d2v  82 

CP - H - 

dt2  ox2 


El 


a2v 

dx2 


=  q(x,t ) 


(5.1) 


where  v(x,t)  is  the  transverse  displacement  of  the  beam,  cp  is  the  mass  density  per 
length,  El  is  the  beam  rigidity,  q(x,t )  is  the  externally  applied  pressure  loading,  t  and 
vindicate  time  and  spatial  axis  along  the  beam  axis.  Rearranging  equation  (5.1)  to 
develop  the  finite  element  formulation  and  the  corresponding  matrix  equations,  the 
averaged  weighted  residual  of  equation  (5.1)  is: 


d2v  d2 

(p - T"  H - T- 

dt 2  dx2 


EI 


d\ 

dx2 


A 


-q 


codx  =  0 


J 


(5.2) 


where  L  is  the  length  of  the  beam  and  co  is  a  test  function.  The  weak  formulation  of 
equation  (5.2)  is  obtained  from  integration  by  parts  twice  for  the  second  term  of  the 
equation.  In  addition,  discretization  of  the  beam  into  a  number  of  finite  elements  gives: 

n  f  -.2  -.2  -.2  ^ 

d  v  ,  r  <3  v  d  co 


f 


+ 


-Veo-M 


V 


da> 

dx 


=  0  (5.3) 
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where  the  slope  is  represented  by: 


Bending  moment: 


Shear  force: 


and  the  Loading: 


0  = 


fdv \ 
UxJ 


M  =  El 


rdV 

Kdx2  y 


V  =  -El 


dx3 


w  = 


f  ^4  \ 
O  V 

dx4 
\ux  J 


(5.4) 


(5.5) 


(5.6) 


(5.7) 


Q‘  represents  the  element  domain  and  n  is  the  number  of  elements  for  the  beam. 


Shape  functions  for  the  spatial  interpolation  of  the  transverse  deflection  v  in 
terms  of  nodal  variables  are  considered.  Elements  are  measured  as  having  two  nodes,  one 
at  each  end  as  shown  in  Figure  17.  The  deformation  of  a  beam  must  have  continuous 
slope  as  well  as  continuous  deflection  at  any  two  neighboring  beam  elements.  To  satisfy 
this  continuity  requirement  each  node  has  both  deflection,  v;  and  slope  0t  as  nodal 


variables.  In  this  case  any  two  neighboring  beam  elements  have  common  deflection  and 
slope  at  the  shared  nodal  point.  This  satisfies  the  continuity  of  both  deflection  and  slope. 
The  Euler-Bemoulli  beam  equation  is  based  on  the  assumption  that  the  plane  normal  to 
the  neutral  axis  before  deformation  remains  normal  to  the  neutral  axis  after  deformation. 


This  assumption  denotes  6  = 
ofx ).  (Kwon  and  Bang,  2000). 


dv 

\dxj 


(i.e.  slope  is  the  first  derivative  of  deflection  in  terms 
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Figure  17.  Two-noded  beam  element.  [After  Kwon  and  Bang,  2000], 


Because  there  are  four  nodal  variables  for  the  beam  element,  we  assume  a  cubic 
polynomial  function  for  v(x) : 

v(x)  =  co  +  qx  +  c2x2  +  c3x3  (5-8) 

The  slope  is  computed  from  equation  (5.8): 

6>(x)  =  c1+2c2x  +  3c3x2  (5-9) 

Evaluation  of  the  deflections  and  slopes  at  the  nodes  yields: 

v(0)  =  c0  =  vl 

m=ct=0t 

2  3  (5-10) 

v(/)  —  Cg  +  Cj/  +  Cjl  +  C3/  —  v2 

0(1)  =  Cj  +  2  c2l  +  3  c3/2  =  02 

where  /  represents  the  element  length.  Solving  equation  (5.10)  for  c.  in  terms  of  the 
nodal  variables  y  and  6i ,  substituting  the  results  into  equation  (5.8)  gives: 


v(x,  t)  =  Efj  (x)  Vj  (t)  +  H2  (x)<9j  (t)  +  H3  (x)  v2  (t)  +  Ha  (x)62  (t)  (5.1 1) 

where  Ht(x)  represents  the  Hermitian  shape  functions,  which  are  of  C1  type,  meaning 

that  both,  v  and  —  are  continuous  between  the  neighboring  elements. 

dx 
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In  general,  C”  type  continuity  means  the  shape  functions  have  continuity  up  to 
the  nth  order  derivative  between  two  neighboring  elements.  The  Hermitian  Shape 
functions  for  the  two  node  beam  element  depicted  in  Figure  17  are: 


3x2  2x3 

*■<*)- 1—+- 


H2  (x)  =  x  - 


2x2 
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3x2  2x3 

w44 
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TT  ,  ,  X  X 

,<x)=“7  F 


(5.12) 


Figure  18.  Hermitian  beam  element.  [From  Kwon  and  Bang,  2000], 


A.  STRAIN  AND  KINETIC  ENERGY  DERIVATION 


Regarding  the  Strain  Energy  due  to  bending,  consider  a  differential  element  of 
beam  of  length  dx ,  where  the  axial  stress  due  to  bending  is: 


Mz 

<7  = - 

/ 


(5.13) 
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where  /  is  the  cross  section  area  of  inertia,  the  Strain  Energy  per  unit  volume  is: 


dw  =  dU  =  —Fdl 
2 


(5.14) 


where  dl  is  the  change  in  length  of  a  fiber  of  the  beam  segment  dx.  The  work  done  by 
the  applied  moment  induces  axial  strain  in  the  fiber,  which  is  stored  as  potential  (strain 
energy). 


Let  F  =  crdA ,  and  dL  =  — L  =  sdx ,  then  the  Srain  Energy  becomes: 

L 


dU  =  ^{<jdA){sdx)  =  ^asdAdx  =  ^asdVol 


(5.15) 


Let  dVol  =  1 ,  then  equation  (5.15)  becomes: 


dU  =  —  crsdVol  =  —  as 
2  2 


(5.16) 


Since  s  = -  and  from  equation  (5.13),  the  Strain  Energy  per  unit  volume 

El 


becomes: 


2 El  I 


(5.17) 


The  total  Strain  Energy  is: 


U=S<IU  =  j  J 

vol  0  A 


1  Mz 


2E\  I 


dAdx 


(5.18) 


Using  equation  (5.5)  and  I  =  ^z2dA,  the  Strain  Energy  stored  in  a  beam  of 


uniform  cross  section  that  has  a  deflected  elastic  curve  v(x)  is: 


El  \( d2v(x,t) 


2  J0{  dx 


(5.19) 
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With  respect  to  the  Kinetic  Energy,  let  the  displacement  for  any  cross  section 
along  the  beam,  to  vibrate  in  a  single  mode  at  its  corresponding  natural  frequency  con : 

v(x,t)  =  Ht{x)  *  Cn  sin (coj  +  0n)  (5.20) 

Where  Cn  %m(ojnt  +  6n )  represents  the  temporal  nodal  coordinates  comprised  of 
displacements  vt(t)  or  slopes  O^t)  as  it  is  implied  from  equation  (5.11).  The  velocity  of 
the  cross  section  becomes: 

v(x,t)  =  Ht(x) *  o onCn  cos {cont  +  0n)  (5.21) 

Let  Cn=  1 ,  then  the  maximum  velocity  is: 

vma  M)  =  Ht{xycon  (5.22) 

Notice  that  vmax  (x)  becomes  a  function  of  xonly.  And  the  Total  Kinetic  Energy 
for  the  beam  is  (mass/length  y ): 

1  L  2 

^max  =  -  WVmax  =  ^\{Hi  W  (5-23) 

Equations  (5.19)  and  (5.23)  were  used  to  derive  the  Strain  and  Kinetic  Energy 
Sensitivities,  respectively  in  the  twenty  element  beam  which  characteristics  are  described 
on  Table  1.  In  order  to  derive  the  Strain  Energies,  in  accordance  with  equation  (5.19),  the 
second  derivatives  of  the  Hermitian  shape  functions  were  multiplied  by  the 
corresponding  displacement  v;  or  rotation  6i  of  the  eigenvector  matrix  (mode  shapes), 

generated  from  the  Build2beams.m  program  in  MATLAB.  Once  the  iterative  process  is 
completed  for  the  number  of  elements,  the  Strain  Energy  matrix  is  assembled  and  it  is  of 
size  mxp,  where  the  m  rows  represent  the  displacement  and  rotations  and  the  p  columns 
represent  the  different  mode  shapes,  as  opposed  to  the  Sensitivity  matrix  where  the  m 
rows  represent  the  mode  shapes  and  the  p  columns  represent  the  different  elements  of  the 
model. 

In  regard  to  the  Kinetic  Energy  calculation,  it  is  implied  from  equation  (5.23),  no 

modifications  were  performed  to  the  Hermitian  shape  functions  which  are  multiplied  by 
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the  natural  frequencies  generated  in  the  Build2beams.m  MATLAB  program.  After  the 
iterative  process  is  completed,  the  Kinetic  Energy  matrix  is  assembled  with  the  same 
characteristics  of  the  Strain  Energy  matrix,  that  is  the  m  rows  representing  the 
displacements  and/or  rotations  and  the  p  columns  representing  the  different  mode  shapes. 

The  following  graphs  show  the  normalized  plots  for  Strain  and  Kinetic  Energies 
distribution,  in  the  same  style  as  the  Sensitivity  distribution  figures  shown  in  the  previous 
chapter. 

Again  each  subplot  corresponds  to  a  specific  mode  shape  and  each  bar  represents 
the  level  of  either  Strain  or  Kinetic  Energy  stored  at  each  element  along  the 
corresponding  mode  shape.  The  red  bars  on  Figure  19  depict  the  Strain  energy 
distribution  along  the  twenty  elements  Base  system  of  the  beam  over  the  first  five  mode 
shapes,  on  this  graph  we  can  see  how  the  Strain  Energy  exactly  correlates  to  the  Stiffness 
Sensitivity  of  the  cantilever  beam  model,  as  it  was  shown  in  Figure  11,  Chapter  IV. 
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The  magenta  color  bars  in  Figure  20  depict  the  Kinetic  Energy  distribution  along 
the  twenty  elements  Base  system  of  the  beam  over  the  first  five  mode  shapes,  on  this 
graph  we  can  conceivably  see  how  the  Kinetic  Energy  is  exactly  correlated  to  the  Mass 
Sensitivity  of  the  cantilever  beam  model,  as  it  was  shown  in  Figure  12,  Chapter  IV. 
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In  the  same  way,  a  comparison  of  the  Strain  and  Kinetic  Energies  to  the  Mass  and 
Stiffness  Sensitivities  is  performed  with  the  ABC’s  applied. 


43 


Figure  21  depicts  the  Strain  Energy  distribution  along  the  twenty  elements  with 
ABC  applied  at  DOF  25  of  the  beam  over  the  first  five  mode  shapes,  in  this  graph  we  can 
see  how  the  Strain  Energy  distribution  is  exactly  correlated  to  the  Stiffness  Sensitivity 
distribution  of  the  cantilever  beam  model  with  the  ABC  applied  at  the  same  DOF,  as  it 
was  showed  in  Figure  13,  Chapter  IV. 
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Figure  22  depicts  the  Kinetic  Energy  distribution  along  the  twenty  elements  with 
ABC  applied  at  DOF  25  of  the  beam  over  the  first  five  mode  shapes,  on  this  graph  we 
can  see  how  the  Kinetic  Energy  distribution  is  exactly  correlated  to  the  Mass  Sensitivity 
of  the  cantilever  beam  model  with  the  ABC  applied  at  the  same  DOF,  as  it  was  showed  in 
Figure  15,  Chapter  IV. 
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It  has  been  proven  for  the  cantilever  beam  modeled,  that  the  Strain  and  Kinetic 
Energy  matrices  are  directly  correlated  to  the  Stiffness  and  Mass  Sensitivity  matrices, 
respectively;  therefore  we  can  conclude  that  the  same  trend  observed  on  the  Stiffness 
Sensitivities  distribution  from  the  last  chapter  applies  for  the  Strain  Energy  distribution, 
as  well  as  the  same  trend  observed  for  the  Mass  sensitivity  distribution  applies  for  the 
Kinetic  Energy  distribution. 

It  is  also  seen  how  the  Stiffness  Sensitivity  matrix  is  directly  correlated  to  the 
flexural  properties  of  the  structure  and  how  the  Mass  Sensitivity  matrix  is  directly 
correlated  to  the  Kinetic  energy  of  the  structure,  which  in  other  words  indicates  it  is  a 
function  of  the  velocity  stored  in  the  structure  at  each  specific  mode  shape.  That  is  why 
the  Figures  depicted  of  Mass  Sensitivity  and  Kinetic  Energy  distributions  showed  low 
sensitivity  at  the  hard  points  of  the  structure  (i.e.,  high  stiffness  areas)  and  higher 
sensitivity  on  the  portions  away  from  the  fixed  points,  where  there  is  minimum  Potential 
energy,  but  maximum  Kinetic  energy. 

This  is  going  to  play  an  important  role  as  it  will  be  seen  in  the  next  chapter,  as  it  is 
a  very  useful  tool  in  the  endeavor  of  predicting  damage  detection  in  the  model. 
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VI.  DAMAGE  DETECTION 


It  was  shown  that  there  exists  a  direct  correlation  between  the  Mass  sensitivity 
with  respect  to  the  Kinetic  Energy  matrices  as  well  as  the  Stiffness  sensitivity  with  regard 
to  the  Strain  Energy  matrices.  Furthermore  with  the  tendencies  observed  and  discussed  at 
the  end  of  the  last  chapter  when  the  ABC’s  were  applied;  several  scenarios  were  run  in 
order  to  find  the  best  choice  of  ABC  coordinates,  such  that  it  could  result  with  the  best 
error  estimate  of  the  structural  model. 

As  the  number  of  elements  is  augmented,  the  finite  element  model  accuracy 
improves.  The  disadvantage  is  that  as  the  number  of  elements  increases,  so  does  the 
computational  time.  The  twenty  element  model  was  used  to  perform  the  error  prediction 
calculation  in  the  present  thesis,  making  use  of  the  first  five  modes  to  perform  the 
analysis,  as  it  was  decided  due  to  the  reliable  convergence,  when  the  number  of  elements 
of  the  model  gets  increased,  as  it  can  be  seen  on  Table  2. 


#  OF  ELEMENTS 

10 

20 

42 

MODE  1 

8.632305 

8.632298 

8.632298 

MODE  2 

54.09948 

54.0978 

54.09769 

MODE  3 

151.5137 

151.4776 

151.4752 

MODE  4 

297.1136 

296.8493 

296.8317 

MODE  5 

491.9194 

490.7656 

490.6868 

MODE  6 

736.9511 

733.2692 

733.0091 

MODE  7 

1033.978 

1024.51 

1023.809 

MODE  8 

1385.246 

1364.733 

1363.099 

MODE  9 

1791.106 

1754.315 

1750.902 

MODE  10 

2226.603 

2193.797 

2187.248 

Table  2.  Frequency  (Hz)  versus  beam  elements. 


In  order  to  run  the  experiments,  two  FEM  beams  were  modeled  in  MATLAB,  one 
being  the  analytical  model  and  the  other  being  the  experimental,  both  models  are  exactly 
the  same  and  correspond  to  the  cantilever  beam  depicted  in  Figure  23,  the  beam 
specifications  are  those  stated  on  Table  1. 
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The  experimental  beam  was  used  to  introduce  a  known  perturbation  at  a  known 
element  as  it  is  requested  from  the  MATLAB  computer  program. 

The  errors  implemented  into  the  experimental  beam  consisted  of  a  10  percent  of 
either  mass  or  stiffness  perturbations. 
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Figure  23.  20  Element  Cantilever  Beam. 


A.  EXPERIMENT 

In  order  to  find  a  set  of  ABC’s,  such  that  we  could  produce  a  set  of  modes  with  an 
improved  sensitivity  distribution,  consequently  a  better  error  estimate;  twenty  four 
different  cases  were  run  to  visualize  the  error  perturbations  introduced  on  the 
experimental  cantilever  beam  model.  Twelve  scenarios  consider  the  analysis  of  one  ABC 
and  one  damaged  element,  and  twelve  scenarios  consider  the  analysis  of  multiple  ABC’s 
and  multiple  damaged  elements.  Table  3  illustrates  the  experiment  set  up  as  follows: 


ONE  ABC/  ONE  DAMAGED 

MULTIPLE  ABC’S/  MULTIPLE 

ELEMENT 

DAMAGED  ELEMENTS 

TYPE  OF  DAMAGE 

NUMBER  OF  CASES 

TYPE  OF  DAMAGE 

NUMBER  OF  CASES 

STIFFNESS 

6 

STIFFNESS 

6 

MASS 

6 

MASS 

6 

Table  3.  Twenty  four  Experiment  Cases. 


On  the  following  cantilever  beam  drawings,  the  pinned  ABC’s  are  represented  by 
the  red  triangles,  the  elements  with  the  stiffness  perturbations  are  represented  in  brown 
color  and  the  elements  with  the  corresponding  mass  perturbations  are  represented  in 
magenta. 
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The  selection  of  the  ABC’s  was  chosen  following  the  sensitivity  patterns 
identified  and  discussed  on  the  last  two  chapters  in  relation  to  the  mass  and  stiffness 
behaviors  with  respect  to  the  Kinetic  and  Strain  Energies  performance,  respectively. 

The  one  ABC/  one  damaged  element  cases  are  analyzed  with  3  different  graphs: 

•  The  cantilever  beam  with  the  graphic  error  representation  of  the  error,  as 
well  as  the  pinned  DOF  where  the  ABC  is  applied, 

•  The  stiffness  or  mass  sensitivity  distribution  along  the  beam  for  the 
specific  case  in  analysis,  and 

•  The  actual  damage  detection  (error  prediction)  of  the  model. 

B.  STIFFNESS  DAMAGE  DETECTION  WITH  ONE  ABC 

1.  Case  1.  Stiffness  Error  at  Element  5,  and  ABC  Pinned  at  DOF  9 


Figure  24.  Stiffness  error  at  element  5,  ABC  pinned  at  DOF  9. 


Figure  25  shows  the  Stiffness  sensitivity  distribution  along  the  beam  represented 
by  the  red  bars,  the  same  figure  also  shows  how  the  sensitivity  is  increased  at  the  DOF’s 
where  the  ABC  is  applied.  Figure  26  represents  the  error  localization  at  those  places 
where  the  blue  bars  are  raised.  The  stem  plot  with  the  red  circle  on  the  top  indicates  the 
magnitude  and  localization  of  the  actual  error. 

Figure  26  indicates,  how  the  damage  detection  is  correctly  identified,  showing 
how  areas  of  strong  sensitivity  are  very  prone  to  detect  the  stiffness  errors  imposed  on  the 
beam.  Additionally  the  increased  sensitivity  due  to  the  pinned  ABC  at  a  DOF  adjacent  to 
the  element  where  the  stiffness  error  perturbation  is  applied  showed  the  effectiveness  in 
the  damage  detection. 
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Figure  25.  Stiffness  sensitivity  distribution  Case  1. 
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Figure  26.  Stiffness  error  prediction  Case  1. 
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2. 


Case  2.  Stiffness  Error  at  Element  5,  and  ABC  Pinned  at  DOF  31 
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Figure  27.  Stiffness  error  at  element  5,  ABC  pinned  at  DOF  3 1 . 


Figure  28  shows  the  Stiffness  sensitivity  distribution  along  the  beam  represented 
by  the  red  bars,  this  figure  shows  how  the  sensitivity  is  increased  at  the  DOF  where  the 
ABC  is  imposed,  which  as  for  this  case  is  located  distant  with  respect  to  the  element 
where  the  stiffness  error  is  applied. 

Even  though  the  sensitivity  bars  in  the  damaged  element  region  are  low,  the  actual 
error  is  located.  This  situation  is  due  to  the  fact  that  the  Base  sensitivity  distribution  (no 
ABC’s  applied)  next  to  the  fixed  end  of  the  beam  is  high,  as  it  is  shown  in  Figure  19,  that 
is  why  the  error  prediction  according  to  Figure  29  still  preserves  the  damage  when  modes 
one  through  three  up  to  modes  one  through  five  are  retained. 

This  case  shows  how  even  when  the  pinned  ABC  is  applied  at  a  DOF  further 
apart  from  the  stiffness  damaged  element,  where  no  increase  in  the  sensitivity  distribution 
is  achieved,  it  still  can  be  able  to  retain  the  error  perturbation,  as  long  as  the  Base 
sensitivity  distribution  is  high  enough  at  that  particular  section  of  the  model. 
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Figure  28.  Stiffness  sensitivity  distribution  Case  2. 
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Figure  29.  Stiffness  error  prediction  Case  2. 
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3. 


Case  3.  Stiffness  Error  at  Element  10,  and  ABC  Pinned  at  DOF  19 


/I 


Figure  30.  Stiffness  error  at  element  10,  ABC  pinned  at  DOF  19. 

Case  3  is  exactly  another  representation  of  case  1,  where  it  can  be  seen  in  Figure 
3 1  the  places  of  high  sensitivity  at  the  pinned  DOF,  which  happens  to  coincide  with  the 
element  among  the  actual  stiffness  error  perturbation. 

This  situation  again  predicts  the  error  in  a  consistent  way  although  not  as  accurate 
as  in  case  1.  When  modes  one  through  three  were  retained,  the  error  was  picked  up 
accurately,  retained  modes  one  through  four  and  one  through  five  bounded  the  error  at 
element  nine,  which  is  still  close  to  the  place  of  the  actual  error  location. 

One  possible  reason  is  due  to  the  fact  that  areas  of  high  strain  energy 
concentration  such  as  the  fixed  end  of  the  cantilever  beam  are  more  likely  to  detect 
stiffness  errors,  thus  the  further  away  the  damage  element  is  from  the  fixed  end  of  the 
cantilever  beam  model  (area  of  high  strain  energy  concentration),  the  more  difficult  the 
damage  detection  becomes.  That  is  why  the  closest  the  damage  element  is  to  the  fixed 
end,  the  better  the  damage  detection  will  work  due  to  the  strain  energy  distribution.  This 
is  fortuitous,  since  the  damage  is  more  likely  to  develop  near  the  fixed  end,  where  the 
maximum  stress  concentration  occurs  (Panday,  1991). 
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Figure  3 1 .  Stiffness  sensitivity  distribution  Case  3. 
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Figure  32.  Stiffness  error  prediction  Case  3. 
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4. 


Case  4.  Stiffness  Error  at  Element  10,  and  ABC  Pinned  at  DOF  31 
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Figure  33 .  Stiffness  error  at  element  1 0,  ABC  pinned  at  DOF  3 1 . 


Case  4  is  another  example  related  to  case  2,  where  Figure  34,  shows  how  the 
sensitivity  is  increased  next  to  the  DOF  where  the  ABC  is  applied,  which  as  for  this  case 
is  located  distantly  with  respect  to  the  element  where  the  stiffness  error  is  exerted, 
although,  this  same  figure  also  shows  how  there  is  still  some  sensitivity  exposure  along 
the  damaged  element. 

The  absence  of  blue  bars  in  Figure  35  over  the  stem  plot,  shows  how  the  error  is 
never  found,  hence  it  confirms  the  notion  that  pinned  ABC  at  a  DOF  further  away  from 
the  element  where  the  stiffness  perturbation  is  applied  is  ineffective,  in  addition,  since  the 
Base  Stiffness  sensitivity  or  Base  Strain  Energy  distributions  are  not  really  high  at 
element  10,  there  are  no  possibilities  that  the  error  can  be  perceived. 

This  case  clearly  demonstrates  how  in  distant  areas  (DOF’s)  where  ABC  is 
applied,  the  sensitivity  is  not  very  likely  to  be  increased;  hence  no  stiffness  error  is 
possible  to  be  located. 
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Figure  34. 

Stiffness  sensitivity  distribution  Case  4. 
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Figure  35.  Stiffness  error  prediction  Case  4. 
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5. 


Case  5.  Stiffness  Error  at  Element  15,  and  ABC  Pinned  at  DOF  31 


Figure  36.  Stiffness  error  at  element  15,  ABC  pinned  at  DOF  3 1 . 


Case  5  is  analogous  to  cases  1  and  3  and  concludes  the  set  of  representations  of 
converged  pinned  ABC  at  the  DOF  location  of  the  stiffness  damaged  element.  This  case 
again  predicts  the  error  in  a  consistent  way,  where  the  error  detection  is  picked  when  the 
second  mode  is  retained  all  the  way  up  to  retained  mode  five  as  it  is  depicted  in  Figure 
38. 

Cases  1,  3,  and  5  prove  the  evidence  found  in  Chapters  IV  and  V  where  in 
particular  for  stiffness  errors,  when  an  ABC  is  particularly  applied  adjacent  to  the 
damaged  element,  the  perturbation  can  be  very  well  detected  and  retained,  due  to  the 
increased  sensitivity  effect  over  the  pinned  ABC  region. 
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Figure  37.  Stiffness  sensitivity  distribution  Case  5. 
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Figure  38.  Stiffness  error  prediction  Case  5. 
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6. 


Case  6.  Stiffness  Error  at  Element  15,  and  ABC  Pinned  at  DOF  11 
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Figure  39.  Stiffness  error  at  element  15,  ABC  pinned  at  DOF  1 1 . 


Case  6  is  the  preceding  and  last  description  of  cases  2  and  4,  where  the  pinned 
ABC  is  applied  at  a  DOF  further  away  from  the  location  of  the  stiffness  damaged 
element. 

Figure  40  shows  how  the  sensitivity  varies  arbitrarily  with  respect  to  the  damage 
element,  hence  is  very  unlikely  to  have  a  reliable  prediction  of  the  error,  as  it  is  the  case 
shown  in  Figure  41,  nevertheless  the  error  happened  to  get  retained  at  modes  one  to  four 
and  modes  one  to  five.  Even  though  this  was  the  situation  aroused  for  this  specific  case, 
there  is  no  consistent  statement  that  this  pinned  ABC  set  up  is  going  to  accurately  predict 
the  stiffness  error  on  the  cantilever  beam  model. 

In  summary,  Cases  2,  4  and  6  ratify  that  when  ABC’s  are  applied  at  a  DOF  further 
away  from  the  element  in  error,  there  is  not  advantage  of  the  increased  stiffness  of  the 
fixed  condition  effect,  therefore  the  sensitivity  is  not  considerably  raised;  thus  no 
stiffness  error  is  likely  to  be  found. 
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Figure  40.  Stiffness  sensitivity  distribution  Case  6. 
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Figure  4 1 .  Stiffness  error  prediction  Case  6. 
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C.  MASS  DAMAGE  DETECTION  WITH  ONE  ABC 

1.  Case  1.  Mass  Error  at  Element  5,  and  ABC  Pinned  at  DOF  11 
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Figure  42.  Mass  error  at  element  5,  ABC  pinned  at  DOF  1 1 . 


Figure  43  shows  the  Mass  sensitivity  distribution  along  the  beam  represented  by 
the  magenta  bars,  this  figure  also  shows  how  the  sensitivity  is  decreased  at  the  DOF 
where  the  ABC  is  applied.  Figure  44  represents  the  error  localization  at  those  places 
where  the  blue  bars  are  raised.  The  stem  plot  with  the  green  circle  on  the  top  indicates  the 
magnitude  and  localization  of  the  actual  mass  error. 

As  for  this  case,  Figure  44  indicates,  how  no  damage  detection  is  identified,  since 
no  blue  bars  are  visible  along  the  twenty  element  beam  model,  showing  how  areas  of  low 
sensitivity  are  not  going  to  detect  the  mass  errors  imposed  on  the  beam. 

Since  the  mass  perturbation  error  is  directly  correlated  to  the  Kinetic  Energy 
distribution  of  the  beam,  the  mass  error  is  a  function  of  the  velocity  (inertia)  of  the  model. 
Consequently  the  possibility  of  detecting  mass  errors  at  DOF’s  where  the  ABC  coincides 
with  the  element  with  the  mass  error  perturbation  are  minimal.  This  artificial 
manipulation  implies  that  the  corresponding  DOF  is  pinned  to  ground,  therefore  no 
displacement  is  very  likely  to  occur. 
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Figure  43.  Mass  sensitivity  distribution  Case  1. 


MODE  1 


0.1 

□ 

-0.1 

0.1 

□ 

-0.1 

0.1 

□ 

-0.1 

0.1 

□ 

-0.1 

0.1 

□ 

-0.1 


1 

I 

u1- 

i 

i  ; 

I 

I 

i 

I 

j 

; 

apt  ; 

1 

! 

i  i 

; 

i 

j 

; 

□ 

2 

4 

6 

8 

10  12 
MODES  1:2 

14 

16 

18 

20 

! 

TV 

i 

i  i 

! 

! 

i 

i 

1 

1 

i 

i  i 

1 

i 

i 

i 

□ 

2 

4 

6 

8 

10  12 
MODES  1:3 

14 

16 

18 

20 

! 

?*! 

i 

i  ; 

i 

! 

i 

i 

i 

1 

i 

i  i 

i 

i 

i 

i 

□ 

2 

4 

6 

.r->. 

8 

10  12 
MODES  1:4 

14 

16 

18 

20 

! 

T*' 

i 

! 

! 

i 

i 

1 

1 

i 

i  i 

1 

i 

i 

i 

□ 

2 

4 

6 

8 

10  12 
MODES  1:5 

14 

16 

18 

20 

1 

i 

i  ; 

! 

i 

! 

i 

i 

■ 

i 

i 

i  i 

i 

1 

i 

i 

□ 

2 

4 

6 

8 

10  12 

14 

16 

18 

20 

Figure  44.  Mass  error  prediction  Case  1. 
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2. 


Case  2.  Mass  Error  at  Element  5,  and  ABC  Pinned  at  DOF  19 
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Figure  45.  Mass  error  at  element  5,  ABC  pinned  at  DOF  19. 


Figure  46  shows  the  mass  sensitivity  distribution  along  the  beam  represented  by 
the  magenta  bars,  in  the  same  figure  it  is  also  seen  how  the  sensitivity  is  pushed  away 
from  the  DOF  where  the  ABC  is  applied,  and  how  the  sensitivity  is  increased  at  those 
DOF  where  the  maximum  displacement  is  present. 

Since  element  five  is  comprised  within  the  area  of  maximum  displacement, 
therefore  good  sensitivity  is  present  at  those  DOF.  The  blue  bars  on  Figure  47  depict  how 
the  error  is  well  retained  in  magnitude  and  location  when  the  first  three  modes  are 
retained.  This  scenario  shows  how  the  mass  perturbation  errors  are  a  function  of  the 
maximum  displacement  along  the  beam. 

Given  that  the  mass  perturbation  error  is  correlated  to  the  Kinetic  energy 
distribution  along  the  beam,  the  areas  where  displacement  is  evident,  will  also  represent 
areas  of  Kinetic  Energy  presence;  therefore  the  error  is  very  likely  to  be  predicted. 
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Figure  46.  Mass  sensitivity  distribution  Case  2. 
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Figure  47.  Mass  error  prediction  Case  2. 
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3. 


Case  3.  Mass  Error  at  Element  10,  and  ABC  Pinned  at  DOF  19 
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Figure  48.  Mass  error  at  element  10,  ABC  pinned  at  DOF  19. 


Case  3  is  analogous  to  case  1 ,  where  the  pinned  DOF  happens  to  overlap  with  the 
element  with  the  actual  mass  error  perturbation.  Figure  49  clearly  shows  how  the 
sensitivity  is  diminished  at  the  DOF  where  the  ABC  is  applied. 

The  absence  of  blue  bars  on  Figure  50  evidently  proves  again  how  no  damage 
detection  is  identified.  This  case  once  again  shows  how  areas  of  low  sensitivity  are  not 
going  to  detect  the  mass  errors  imposed  on  the  beam.  Again,  the  lack  of  displacement 
effect  due  to  the  ABC  pinned  at  a  DOF  adjacent  to  the  element  mass  perturbation 
demonstrates  how  the  mass  error  prediction  is  not  likely  to  be  successful. 


65 


Mode  1 


- 

- 1 - 

L 

- 1 - 

i 

- 1 - 

L 

- 1 - 

'JL'-.  — 

- 1 - 

— i - 

- 1 - 1 - 1 — 

. . 

M 

n 

□ 

2 

6 

s 

10 

Mode  2 

12 

14 

16 

16 

20 

- 

, n  n  f 

i  i 

1  1  . 

i 

.  n  i  n  ■  „  ' 

I] 

□ 

□ 

2 

z 

[ 

6 

s 

10 

Mode  3 

12 

14 

16 

10 

20 

- 

■  1 

rrri 

1 

1  ■ . 

i 

-  1 

iii.. 

„  l 

r 

□ 

2 

z 

[ 

6 

S 

10 

Mode  4 

12 

14 

16 

10 

20 

- 

. .  i 

1 1 

nrm 

□ 

i 

f  n 

i  i 

i . . 

i 

.  n  f 

uj 

□ 

□ 

2 

z 

[ 

6 

0 

10 

Mode  5 

12 

14 

16 

10 

20 

- 

.  1 

1 1 

i . ;  i 

□ 

u 

i 

.  1 

D 

.  i 

1  n  .  [ 

L 

□ 

2 

z 

[ 

6 

0 

10 

12 

14 

16 

10 

20 

Figure  49.  Mass  sensitivity  distribution  Case  3. 
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Figure  50.  Mass  error  Prediction  Case  3. 
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4. 


Case  4.  Mass  Error  at  Element  10,  and  ABC  Pinned  at  DOF  31 


Figure  5 1 .  Mass  error  at  element  10,  ABC  pinned  at  DOF  3 1  and  DOF  13. 


Case  4  happened  to  be  a  very  sturdy  scenario  on  finding  the  mass  damage 
detection.  It  was  assumed  that  the  pinned  ABC  that  would  generate  one  of  the  highest 
mass  sensitivities  measurements  was  that  at  DOF  31.  Figure  52  shows  the  mass 
sensitivity  distribution  along  the  beam  model  for  the  mentioned  pinned  DOF.  But  Figure 
53  showed  how  the  damage  detection  was  never  found,  hence  forcing  to  run  the  scenario 
for  a  variety  of  pinned  ABC  conditions. 

After  several  runs  pinned  ABC  at  either  DOF  13  or  DOF  15  (dark  red  triangle  in 
Figure  51),  proved  to  predict  the  mass  error  perturbation  exactly  in  the  same  approach  as 
it  is  depicted  on  Figure  55,  where  the  magnitude  and  location  of  the  error  is  found  when 
the  first  three  modes  are  retained  and  hold  the  error  up  to  the  five  retained  modes. 

DOF  13  happens  to  stay  literally  to  the  same  distance  with  respect  to  element  10, 
as  DOF  3 1  is,  but  the  error  prediction  pattern  was  totally  different  when  analyzed,  where 
DOF  13  proved  to  do  an  acceptable  error  prediction,  but  the  pinned  ABC  at  DOF  31  did 
not  work  as  efficiently  as  was  expected. 

Although  the  evidence  found  in  Chapters  IV  and  V  still  holds  true,  that  is,  areas  of 
higher  displacements  are  very  likely  to  predict  better  the  mass  errors,  it  is  cumbersome, 
how  it  will  not  work  for  every  case  where  the  higher  displacements  are  evident,  making 
this  situation  a  problem  of  major  concern  where  mass  error  detection  is  of  interest. 
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Figure  52.  Mass  sensitivity  distribution  Case  4  at  DOF  31. 
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Figure  53.  Mass  error  prediction  case  4  at  DOF  3 1 . 
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Figure  54.  Mass  sensitivity  Case  4  at  DOF  13. 
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Figure  55.  Mass  error  prediction  case  4  at  DOF  13. 
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5. 


Case  5.  Mass  Error  at  Element  15,  and  ABC  Pinned  at  DOF  31 


Figure  56.  Mass  error  at  element  15,  ABC  pinned  at  DOF  31. 


Case  5  is  the  last  representation  of  cases  1  and  3,  where  the  pinned  DOF  lies 
adjacent  to  the  element  with  the  actual  mass  error  perturbation.  Figure  49  clearly  shows 
how  the  sensitivity  is  diminished  at  the  DOF  where  the  ABC  is  applied. 

Figure  50  proves  again  how  no  damage  detection  is  identified.  This  case 
consolidates  the  reflection  of  how  areas  of  low  sensitivity  are  not  going  to  detect  the  mass 
errors  imposed  on  the  beam.  Again,  the  lack  of  displacement  effect  due  to  the  ABC 
pinned  at  a  DOF  adjacent  to  the  element  mass  perturbation  demonstrates  how  the  mass 
error  prediction  is  never  going  to  be  successful. 
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Figure  57. 

Mass  sensitivity  distribution  Case  5. 
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Figure  58.  Mass  error  prediction  Case  5. 
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6. 


Case  6.  Mass  Error  at  Element  15,  and  ABC  Pinned  at  DOF  23 


Figure  59.  Mass  error  at  element  15,  ABC  pinned  at  DOF  23. 


Case  6  represented  a  very  similar  situation  to  the  one  experienced  in  case  4. 
Initially  DOF  19  was  selected  as  the  pinned  ABC  condition  to  analyze,  as  it  was  assumed 
to  generate  a  good  Mass  sensitivity  distribution,  but  when  the  experiment  was  evaluated, 
it  did  not  find  the  error  accurately.  Hence  as  in  case  4,  this  case  was  run  for  a  variety  of 
pinned  ABC’s  conditions. 

The  analysis  was  organized  in  two  parts.  The  first  with  the  ABC  pinned  to  the  left 
hand  side  of  the  mass  damaged  element,  giving  the  boundary  condition  effect  of  a 
cantilever  beam  model.  After  several  runs  pinned  ABC  at  DOF  23,  was  the  best  ABC 
configuration  able  to  find  the  mass  perturbation  error,  as  it  is  shown  on  Figure  61,  where 
modes  one  to  three  and  one  to  four,  retained  the  error  in  magnitude  and  location,  mode 
five,  lost  the  actual  error  detection,  but  still  found  some  mass  perturbation  at  elements 
fourteen  and  sixteen  of  the  experimental  cantilever  beam. 

The  second  part  of  the  analysis  was  with  the  ABC  set  up  to  the  right  hand  side  of 
the  mass  damaged  element,  giving  the  boundary  condition  effect  of  a  clamped-  clamped 
beam.  The  best  error  prediction  was  found  at  pinned  ABC  at  DOF  35  and  37,  which 
found  the  error  at  mode  five  only,  the  rest  of  the  DOF  did  not  find  the  error  at  all,  and  no 
error  prediction  pattern  proved  to  exist. 
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Figure  60.  Mass  sensitivity  distribution  Case  6. 
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Figure  6 1 .  Mass  error  prediction  Case  6. 
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D.  STIFFNESS  DAMAGE  DETECTION  WITH  MULTIPLE  ABC’S 

With  the  usefulness  of  the  ABC’s  proven  for  a  single  pinned  ABC  and  one  single 
damaged  element,  and  using  as  guidance  the  different  patterns  found  according  to  the 
different  configuration  of  the  ABC’s;  it  is  inferred  that  by  the  use  of  more  than  one  ABC, 
the  system  is  going  to  become  less  underdetermined,  therefore  it  is  going  to  be  able  to 
better  predict  the  error  perturbations  imposed  in  the  cantilever  beam  model. 

The  main  goal  is  to  configure  the  minimum  number  of  ABC’s  such  that  high 
enough  sensitivity  levels  are  attained  throughout  the  structure,  hence  being  able  to  have 
reliable  damage  detection  capability  at  all  times  and  at  any  beam  location,  considering 
that  in  a  real  life  structure  the  location  of  the  damaged  elements  are  very  likely  to  be 
unknown.  This  goal  is  going  to  be  explored  by  the  use  of  three  and  five  ABC’s 
configurations  along  the  cantilever  beam,  for  one,  two,  and  four  damaged  elements  in 
order  to  test  the  ability  of  predicting  the  different  stiffness  errors  in  a  reliable  form. 

The  following  Figures  are  intended  to  represent  damage  detection  for  a  different 
number  and  locations  of  elements  with  an  imposed  stiffness  error,  by  the  use  of  multiple 
ABC’s  along  the  modeled  structure.  The  sensitivity  distribution  graphs  for  the  following 
scenarios  were  omitted,  since  it  is  already  understood  that  the  sensitivity  levels  are  going 
to  get  increased  a  those  DOF’s  where  the  pinned  ABC’s  are  applied. 

The  pinned  ABC’s  are  represented  by  the  red  triangles,  the  elements  with  the 
stiffness  perturbations  are  represented  in  brown  color.  The  stem  plots  with  the  blue 
circles  on  the  top  indicate  the  magnitude  and  localization  of  the  actual  stiffness  errors. 

1.  One  Stiffness  Damaged  Element  with  Three  ABC’s 

Figure  62  shows  the  cantilever  beam  representation  for  10  percent  stiffness 
damage  at  element  ten,  with  a  three  ABC’s  set  up  applied  at  DOF’s  15,  29,  and  41 : 
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Figure  62.  Stiffness  error  at  element  10,  ABC’s  pinned  at  DOF’s  15,  29,  41. 


The  error  prediction  for  the  above  configuration  is  the  following: 
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Figure  63.  Stiffness  error  prediction  at  element  10  with  3  ABC’s. 


Figure  64  shows  the  cantilever  beam  representation  for  a  10  percent  stiffness 
damage  at  element  five,  and  with  the  same  three  ABC’s  set  up: 
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Figure  64.  Stiffness  error  at  element  5.  ABC’s  pinned  at  DOF’s  15,  29,  41. 


The  error  prediction  for  the  above  configuration  is  the  following: 
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Figure  65.  Stiffness  error  prediction  at  element  5  with  3  ABC’s. 


2.  Two  Stiffness  Damaged  Elements  with  Three  ABC’s 

Figure  66  depicts  the  same  three  ABC  configuration,  but  now  for  a  10  percent 
stiffness  damaged  at  elements  5  and  18: 
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Figure  66.  Stiffness  error  at  elements  5,  and  18.  ABC’s  pinned  at  DOF’s  15,  29,  41. 


The  error  prediction  for  the  above  configuration  is  the  following: 


Figure  67.  Stiffness  error  prediction  at  elements  5,  and  18.  3  ABC’s. 


3.  Four  Stiffness  Damaged  Elements  with  Three  ABC’s 

Figure  68  depicts  the  standard  three  ABC’s  set  up  chosen  for  these  experiments, 
for  a  10  percent  stiffness  damage  at  elements  5,  10,  15,  and  20: 
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Figure  68.  Stiffness  error  at  elements  5,  10,  15,  and  20.  ABC’s  pinned  at  DOF’s  15, 

29,41. 


The  error  prediction  for  the  above  configuration  is  the  following: 


Figure  69.  Stiffness  error  prediction  at  elements  5,  10,  15,  and  20.  3  ABC’s. 


4.  One  Stiffness  Damaged  Element  with  Five  ABC’s 

Figure  70  shows  the  cantilever  beam  representation  for  10  percent  stiffness 
damage  at  element  ten,  with  a  five  ABC’s  set  up  applied  at  DOF’s  7,  15,  23,  33  and  41: 
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Figure  70.  Stiffness  error  at  element  10.  ABC’s  pinned  at  DOF’s  7,  15,  23,  33,  41. 

The  error  prediction  for  the  above  configuration  is  the  following: 
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Figure  71.  Stiffness  error  prediction  at  element  10.  5  ABC’s. 


Figure  72  shows  the  cantilever  beam  representation  for  a  10  percent  stiffness 
damage  at  element  five,  and  with  the  same  five  ABC’s  configuration: 
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Figure  72.  Stiffness  error  at  element  5.  ABC’s  pinned  at  DOF’s  7,  15,  23,  33,  41. 


The  error  prediction  for  the  above  configuration  is  the  following: 
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Figure  73.  Stiffness  error  prediction  at  element  5.  5  ABC’s. 


5.  Two  Stiffness  Damaged  Elements  with  Five  ABC’s 

Figure  74  depicts  the  five  ABC’s  set  up,  for  a  10  percent  stiffness  damaged  at 
elements  5  and  18: 
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Figure  74.  Stiffness  error  at  elements  5  and  18.  ABC’s  pinned  at  DOF’s  7,  15,  23,  33, 

41. 

The  error  prediction  for  the  above  configuration  is  the  following: 


Figure  75.  Stiffness  error  prediction  at  elements  5  and  18.  5  ABC’s. 


6.  Four  Stiffness  Damaged  Elements  with  Five  ABC’s 

Figure  76  depicts  the  five  ABC’s  set  up,  but  now  for  a  10  percent  stiffness 
damage  at  elements  5,  10,  15,  and  20: 
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Figure  76.  Stiffness  error  at  elements  5,  10,  15,  and  20.  ABC’s  pinned  at  DOF’s  7, 

15,23,33,41. 

The  error  prediction  for  the  above  configuration  is  the  following: 
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Figure  77.  Stiffness  error  prediction  at  elements  5,  10,  15,  and  20.  5  ABC’s. 
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E. 


MASS  DAMAGE  DETECTION  WITH  MULTIPLE  ABC’S 


The  following  six  scenarios  represent  the  damage  detection  capabilities  of 
multiple  ABC’s  for  one,  two,  and  four  mass  damaged  elements.  In  the  following  Figures 
the  pinned  ABC’s  are  represented  by  the  red  triangles,  the  elements  with  the  mass 
perturbations  are  represented  in  magenta  color.  The  stem  plots  with  the  yellow  circles  on 
the  top  indicate  the  magnitude  and  localization  of  the  actual  mass  errors. 


1.  One  Mass  Damaged  Element  with  Three  ABC’s 


Figure  78  shows  the  cantilever  beam  representation  for  10  percent  mass  damage 
at  element  eight,  with  a  three  ABC’s  configuration  applied  at  DOF’s  15,  29,  and  41: 


Figure  78.  Mass  error  at  element  8,  ABC’s  pinned  at  DOF’s  15,  29,  41. 


The  error  prediction  for  the  above  configuration  is  the  following: 
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2. 


Two  Mass  Damaged  Elements  with  Three  ABC’s 


Figure  80  shows  the  cantilever  beam  representation  for  10  percent  mass  damage 
at  elements  five  and  fifteen,  with  a  three  ABC’s  set  up  applied  at  DOF’s  15,  29,  and  41 : 
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Figure  80.  Mass  errors  at  elements  5  and  15,  ABC’s  pinned  at  DOF’s  15,  29,  41. 

The  error  prediction  for  the  above  configuration  is  the  following: 


Figure  81.  Mass  errors  prediction  at  elements  5  and  15.  3  ABC’s. 
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Figure  82  shows  the  cantilever  beam  representation  for  10  percent  mass  damage 
at  elements  ten  and  eighteen,  with  a  three  ABC’s  set  up  applied  at  DOF’s  15,  29,  and  41: 


Figure  82.  Mass  errors  at  elements  10  and  18,  ABC’s  pinned  at  DOF’s  15,  29,  41. 


The  error  prediction  for  the  above  configuration  is  the  following: 


Figure  83.  Mass  errors  prediction  at  elements  10  and  18.  3  ABC’s. 
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3. 


Four  Mass  Damaged  Elements  with  Three  ABC’s 


Figure  84  shows  the  cantilever  beam  representation  for  10  percent  mass  damage 
at  elements  five,  ten,  fifteen,  and  twenty,  with  a  three  ABC’s  configuration  applied  at 
DOF’s  15,  29,  and  41: 
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Figure  84.  Mass  errors  at  elements  5,  10,  15,  and  20.  ABC’s  pinned  at  DOF’s  15,  29, 

41. 


The  error  prediction  for  the  above  configuration  is  the  following: 
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Figure  85.  Mass  errors  prediction  at  elements  5,  10,  15,  and  20.  3  ABC’s. 
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4. 


One  Mass  Damaged  Element  with  Five  ABC’s 


Figure  86  shows  the  cantilever  beam  representation  for  10  percent  mass  damage 
at  element  eight,  with  a  five  ABC’s  configuration  applied  at  DOF’s  7,  15,  23,  33,  and  41: 
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Figure  86.  Mass  error  at  element  8.  ABC’s  pinned  at  DOF’s  7,  15,  23,  33,  and  41. 


The  error  prediction  for  the  above  configuration  is  the  following: 


Figure  87.  Mass  error  prediction  at  element  8.  5  ABC’s. 
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5. 


Two  Mass  Damaged  Elements  with  Five  ABC’s 


Figure  88  shows  the  cantilever  beam  representation  for  10  percent  mass  damage 
at  elements  seven  and  eleven,  with  a  five  ABC  configuration  applied  at  DOF’s  7,  15,  23, 
33,  and  41: 


Figure  88.  Mass  error  at  elements  7  and  1 1  ABC’s  pinned  at  DOF’s  7,  15,  23,  33,  and 

41. 


The  error  prediction  for  the  above  configuration  is  the  following: 


Figure  89.  Mass  error  prediction  at  elements  7  and  11.5  ABC’s. 
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Figure  90  shows  the  cantilever  beam  representation  for  10  percent  mass  damage 
at  elements  five  and  nine,  with  a  five  ABC’s  set  up  applied  at  DOF’s  7,  15,  23,  33,  and 
41: 
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Figure  90.  Mass  error  at  elements  5  and  9.  ABC’s  pinned  at  DOF’s  7,  15,  23,  33,  and 

41. 


The  error  prediction  for  the  above  configuration  is  the  following: 


Figure  91.  Mass  error  prediction  at  elements  5  and  9.  5  ABC’s. 
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6. 


Four  Mass  Damaged  Elements  with  Five  ABC’s 


Figure  92  shows  the  cantilever  beam  representation  for  10  percent  mass  damage 
at  elements  three,  seven,  eleven,  and  seventeen,  with  a  five  ABC’s  set  up  applied  at 
DOF’s  7,  15,  23,  33,  and  41: 


Figure  92.  Mass  error  at  elements  3,  7,  11,  and  17.  ABC’s  pinned  at  DOF’s  7,  15,  23, 

33,  and  41. 

The  error  prediction  for  the  above  configuration  is  the  following: 


Figure  93.  Mass  error  prediction  at  elements  3,  7,  11,  and  17.  5  ABC’s. 
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Figure  94  shows  the  cantilever  beam  representation  for  10  percent  mass  damage 
at  elements  five,  nine,  fifteen  and  eighteen,  with  a  five  ABC’s  set  up  applied  at  DOF’s  7, 
15,  23,  33,  and  41: 


Figure  94.  Mass  error  at  elements5,  9,  15,  and  18.  ABC’s  pinned  at  DOF’s  7,  15,  23, 

33,  and  41. 


The  error  prediction  for  the  above  configuration  is  the  following: 
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Figure  95.  Mass  error  prediction  at  elements  5,  9,  15,  and  18.  5  ABC’s. 
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F.  DISCUSSION  OF  RESULTS 

The  sensitivity  of  the  natural  frequencies  of  the  structure  to  changes  in  stiffness 
and  mass  were  calculated  from  the  finite  element  analysis. 

The  twelve  cases  analyzed  for  the  changes  in  stiffness;  six  for  one  ABC  and  one 
damaged  element,  and  six  for  multiple  ABC’s  and  multiple  damaged  number  of  elements, 
demonstrated  a  very  specific  behavior. 

For  the  one  ABC  and  one  damaged  element  cases,  the  highest  probability  of 
detecting  the  damage  proved  to  be  when  the  stiffness  error  was  the  closest  to  the  fixed 
end  of  the  cantilever  beam,  which  happens  to  be  the  highest  exposed  part  to  stiffness 
failure  due  to  the  increased  stress  concentration  in  that  area.  The  further  the  stiffness  error 
perturbation  was  translated  to  the  free  end,  the  less  likely  it  was  detected,  hence  the  use  of 
the  pinned  ABC  at  a  DOF  adjacent  to  the  element  in  error,  highly  increased  the 
possibilities  of  detecting  the  damage. 

When  the  pinned  ABC  was  applied  to  a  DOF  further  away  with  respect  to  the 
location  of  the  damage  element,  the  error  was  never  found,  unless  it  was  located  very 
close  to  the  fixed  end  of  the  cantilever  beam,  where  the  Stiffness  Base  sensitivity  assisted 
in  the  detection  of  it. 

On  the  other  hand,  the  multiple  ABC  and  multiple  damaged  number  of  elements 
cases,  significantly  improved  the  damage  detection  capacity.  For  the  3  ABC’s 
configuration,  the  damage  was  successfully  found  up  to  two  damaged  elements,  when 
four  damaged  elements  were  introduced,  the  system  accurately  predicted  three  out  of  the 
four  damaged  elements,  and  the  fourth  element  in  error  although  it  was  not  detected,  was 
very  well  bounded,  as  it  is  depicted  in  Figure  69.  The  five  ABC  configuration 
successfully  found  the  damage  imposed  over  the  elements  in  error.  Further  experiments 
(not  depicted  in  the  present  thesis)  with  more  than  only  four  damaged  elements 
successfully  retained  the  errors  imposed  over  the  system  structure;  this  explains  how  due 
to  the  stable  behavior  of  the  ABC’s  when  to  stiffness  errors  is  concerned  and  making  the 
problem  less  underdetermined  with  the  incorporation  of  more  ABC’s,  the  error  is  very 
likely  to  be  found. 
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The  twelve  cases  analyzed  for  the  changes  in  mass;  six  for  one  ABC  and  one 
damaged  element,  and  six  for  multiple  ABC’s  and  multiple  damaged  number  of  elements, 
did  not  prove  the  same  consistencies  as  the  ones  found  in  the  stiffness  errors. 

Since  the  mass  errors  are  related  to  the  inertia  of  the  system,  the  pinned  ABC’s 
were  set  up  in  such  an  approach  that  were  able  to  produce  a  considerable  displacement  at 
the  location  of  the  damaged  element,  nevertheless,  while  this  reflection  holds  true,  not  all 
of  the  different  pinned  ABC  combinations  at  those  DOF’  that  produced  higher  kinetic 
energy  seemed  to  detect  the  mass  error  perturbation  successfully. 

For  the  one  ABC  and  one  damaged  element  cases,  when  the  pinned  ABC  was  set 
up  to  emulate  the  beam  model  as  if  it  were  a  clamped-clamped  beam,  the  error  was 
sparsely  detected  and  with  no  pattern  defined;  but  the  best  error  predictions  were  found 
when  the  ABC  was  configured  to  emulate  the  beam  model  as  a  cantilever  beam,  exerting 
the  maximum  inertia  (displacement)  at  the  damaged  element. 

On  regard  to  the  multiple  ABC’s  and  multiple  damaged  number  of  elements 
cases,  the  pattern  found  for  the  one  ABC  and  one  mass  damage  cases,  stills  holds  the 
same.  It  was  discovered  that  even  when  the  configuration  of  more  ABC  helped  to  reduce 
the  undetermined  nature  of  the  finite  element  problem,  therefore  being  able  to  predict  a 
few  more  elements  in  error,  no  specific  pattern  was  found.  Nevertheless  it  was 
recognized  that  the  mass  errors  are  a  very  strong  function  of  the  kinetic  energy  of  the 
system.  Figures  81,  85,  87,  89,  and  93,  show  how  the  errors  located  at  low  kinetic  energy 
regions  were  not  detected,  whereas,  Figures  83,  91,  and  95  show  how  the  errors  were 
accurately  predicted  when  a  smart  combination  of  ABC  produced  areas  of  high  kinetic 
energy  at  the  regions  where  the  damaged  elements  were  located. 

On  the  following  chapter  an  Optimization  approach  is  develop  in  order  to  find  the 
mass  error  in  a  cantilever  beam.  As  it  will  be  seen,  the  MATLAB  built  in  optimization 
function,  proves  to  work  as  a  reliable  technique. 
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VII.  OPTIMIZATION  APPROACH 


The  purpose  of  numerical  optimization  is  to  aid  in  rationally  searching  for  the 
design  variables  to  satisfy  constraints  (Vanderplaats,  1984).  That  is,  to  match  the 
dynamic  response  of  a  finite  element  model  to  that  of  the  experimental  response  of  the 
system  of  interest. 

An  optimization  program  was  developed  to  predict  the  mass  error 
detection/prediction  on  the  cantilever  beam  model. 

A.  OPTIMIZATION  THEORY 

The  general  problem  statement  for  a  non-linear  constrained  optimization  problem 
is: 

To  minimize  f(x )  Objective  Function 

Subject  to: 

Sj(x)  -  0  (7.1) 

where  j  =  \,m  Inequality  Constraints 

K  =  0  (7.2) 

where  k  =  1,  p  Equality  Constraints 

x;<x,.<x“  (7.3) 


where  i  =  \,n  Side  Constraints 


The  design  variables  vector  is  represented: 
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\ 

x2 


(7.4) 


x„ 


The  constraint  equations  are  used  to  bound  possible  solution  combinations  to 
meet  the  requirements  of  the  given  situation.  Most  optimization  algorithms  require  that 
an  initial  set  of  design  variables,  x°  be  specified.  Beginning  from  this  starting  point,  the 
design  is  updated  iteratively.  This  process  can  be  written  as: 

xq=xq~1+asq  (7.5) 

where  q  represents  the  iteration  number  and  5  is  a  vector  search  direction  in  the  design 
space.  The  scalar  quantity  a  defines  the  distance  that  we  wish  to  move  in  direction  5. 
The  updated  values  of  x  are  used  to  calculate  the  value  of  the  objective  function  for  each 
iteration.  The  search  direction  is  chosen  to  decrease  the  value  of  the  objective  function 
while  staying  within  the  constraints  of  the  system.  The  process  continues  until  the 
objective  function  converges  to  a  local  minimum  and  the  optimal  solution  is  localized. 

1.  Modal  Assurance  Criterion  (MAC) 

The  modal  assurance  criterion  is  a  scalar  constant  relating  the  causal  relationship 
between  two  modal  vectors.  The  MAC  is  used  as  a  means  to  compare  analytically  and 
experimentally  obtained  vibrational  mode  shapes.  The  MAC  will  have  a  value  between  0 
and  1 .  A  value  of  0  indicates  that  the  two  modal  vectors  are  not  correlated  while  a  value 
of  1  indicates  the  modal  vectors  are  correlated.  The  method  of  calculating  the  MAC  for 
comparing  the  ith  mode  of  the  analytical  system  to  the  ith  model  of  the  experimental 
system  is: 
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MAC  = 


(7.6) 


B.  PROBLEM  FORMULATION 

The  optimization  problem  was  analyzed  according  to  the  cantilever  beam 
displayed  on  Figure  96,  with  the  same  beam  specifications  described  on  Table  1  from 
Chapter  IV.  The  beam  consisted  of  42  elements,  each  1  inch  in  length  and  DOF  1  was  at 
the  fixed  end  of  the  cantilever  beam;  this  corresponded  to  FE  model  element  quantity  and 
length.  A  Series  336  FLEXCEL  ICP  accelerometer  (serial  number  10860)  was  threaded 
into  position  at  node  42  of  the  beam  and  wired  into  Channel  2  on  a  DACTRON  Focus 
front  end  digital  signal  processor  (DSP).  An  excitation  was  applied  by  a  PCB  Series  086 
B03  impact  hammer  (serial  number  269),  which  was  wired  into  Channel  1  on  the  DSP. 
The  accelerometer  and  force  hammer  were  calibrated  using  a  pendulous  test;  the 
sensitivity  adjustment  was  applied  in  the  set-up  of  RT  Pro  Focus  5.57  software. 
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Figure  96.  Experimental  beam  set  up. 
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The  experiment  was  performed  in  two  segments:  The  analytical  and  the 
experimental.  In  the  analytical,  the  Build2 beams. m  MATLAB  code  was  used  to  generate 
the  original  finite  element  beam  model  set  of  original  eigenvalues  and  eigenvectors,  then 
the  damage  (mass  perturbation)  used  during  the  real  experiment,  was  simulated  in 
MATLAB,  generating  the  modified  set  of  eigenvalues  and  eigenvectors;  subsequently 
these  modal  parameters  were  evaluated  into  the  MATLAB  optimization  program.  The 
mass  error  perturbation  was  imposed  on  element  35  of  the  cantilever  beam,  with  a  10% 
mass  excess  with  respect  to  the  total  beam  mass.  The  finite  element  model  was  analyzed 
from  modes  one  to  eight. 

The  second  (experimental)  segment  consisted  of  the  actual  measured  data  from 
the  laboratory.  It  was  created  after  tapping  the  beam  shown  in  figure  62,  producing  the 
corresponding  Frequencies  Response  Forces  (FRF)  by  means  of  the  softwares: 
DACTRON  RT  Pro  Focus  5.57,  and  ME’scopeVES;  node  42  remained  the  reference  as 
the  load  cell  roved  from  one  node  to  the  next  allowing  for  the  measurement  of  the 
response  at  each  node. 

The  modal  data  extracted  out  of  the  ME’scopeVES  software  comprises  the  natural 
frequencies,  damping  ratios  and  magnitudes  and  phases  of  the  different  mode  shapes  of 
the  described  cantilever  beam  with  the  mass  perturbation  added. 

1.  The  RT  Pro  Focus  5.57  software  “Real-time”  was  configured  to  measure  3200 
spectral  lines,  8192  points,  with  a  delta  T  of  166.7ps  over  the  frequency  range  0-2400Hz. 
The  frequency  range  of  0-2400  Hz  was  chosen  because  it  covered  the  first  10  modes  of 
system  and  signal  resolution  was  sufficient  for  data  acquisition.  The  excitation  signal 
proved  to  be  clean  and  thus  no  window  was  used  for  data  measuring. 


2.  Channel  set-up 

Channel  1  (Excitation)  Channel  2(Response) 
Max  Volts  (mV):  0.1  0.3 

Quantity:  Force  Accel. 
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EU:  lbf  gn 

mv/EU:  8.67  104.383 

Coupling:  ICP  AC  0.7  Hz  ICPAC0.7Hz 

Sensitivity  Adjustment:  0  0 

3.  Trigger  set-up 
Source:  Analog  input 

Run  Mode:  Manual  Arm  every  frame 

Input:  Channel  1 

Slope:  Bi-polar 

Level  (%):  1,  Level  (V):  0 

Pre/Post  Points  (-/+):  -10 

Pre/Post  Time  (-/+):  -1.67ps 

4.  Average  set-up: 

Type:  Linear 
Domain:  Frequency 

Frames:  3  (Each  node  was  excited  3  times  and  an  average  taken  and  saved.) 
Accept/Reject:  Manual  Accept/Reject  every  frame. 

(The  user  rejects  double  taps,  under  powered  or  overloaded  signals.) 

5.  Modal  Coordinate  set-up: 

Auto  increment:  ON 
Rove:  Excitation 

Point  increment:  1 
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Export:  UFF  text  format,  frequency  response. 


Once  the  experimental  modal  parameters  were  extracted,  a  Complex  to  Real 
Conversion  was  performed  on  the  mode  shapes.  Complex  ones,  as  derived  from  analysis 
of  test  data,  and  real  ones,  such  as  are  traditionally  produced  from  theoretical  analyses  in 
the  absence  of  a  detailed  knowledge  of  the  nature,  extent  and  distribution  of  damping 
(Ewins,  2000).  The  corresponding  information  was  exported  to  MATFAB  as  a  “.txt  file” 
such  that  the  corresponding  experimental  eigenvalues  and  eigenvectors  could  be  read  in  a 
matrix  form. 

The  modal  to  spatial  transformation  was  obtained  by  the  use  of  the  “Simple 
Method”  (Ewins,  2000);  through  this  method  the  modulus  and  the  phase  of  the  modal 
eigenvectors  were  obtained  from  the  MEscopesVES  software.  The  modulus  being 
interpreted  as  the  square  root  of  the  sum  of  the  squares  of  the  real  and  imaginary  parts  of 
the  modal  analysis  data,  and  the  phase  was  assigned  to  each  node  either  as  0  or  180 
degrees,  which  represented  the  positive  or  negative  value  of  the  displacement  part  of  the 
experimental  eigenvector  matrix.  This  process  was  performed  within  the  MEscopeVES 
software  once  the  shapes  table  was  generated,  by  selecting  “Complexity  plot”  from  the 
“Display”  drop  down  menu,  obtaining  a  ‘Starbust  plot’  -  Argand  Diagram  (Ewins,  2000), 
then  the  modulus  and  phase  table  was  automatically  obtained  after  clicking  on  the 
“Normalize  Data”  button,  now  being  ready  to  import  the  table  as  the  “.txt  file”  into 
matlab.  Figure  97  depicts  the  Starbust  plot  obtained  from  the  10%  mass  error  exerted  at 
element  35  of  the  cantilever  beam. 
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Figure  97.  Starbust  Plot-  Argand  Diagram. 


C.  OPTIMIZATION  APPROACH  VIA  MATLAB 

The  actual  optimization  program  was  run  with  the  use  of  a  built  in  MATLAB 
function  named  “fmincon”  which  permitted  the  finding  of  the  error  perturbation  on  the 
beam,  “fmincon”  aids  in  finding  a  constrained  minimum  of  a  given  Objective  Function 
of  several  design  variables  starting  at  an  initial  estimate  (xo).  The  fmincon  algorithm 
employs  the  SQP  implementation,  which  consists  of  three  stages.  The  first  stage  updates 
the  Hessian  matrix  via  the  quasi-  Newton  method,  where  the  Hessian  matrix  is  calculated 
using  the  Broydon-  Fletcher  Goldgarb-Shanno  (BFGS)  method,  which  is  the  most  used 
with  inexact  line  searches.  The  second  stage  performs  the  Quadratic  Programming 
Solution,  used  to  calculate  the  parameter  changes,  by  the  calculation  of  a  feasible  point 
and  then  the  generation  of  an  iterative  sequence  of  feasible  points  that  converge  to  the 
solution.  The  last  stage  essentially  performs  the  line  search  and  generates  the  merit 
function.  The  “fmincon”  MATLAB  set  up  used  was:  x  =  fmincon 
(fun,xo,A,  b,Aeq,  beq,  lb,  ub) 
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The  “fmincon”  needs  as  input  arguments  the  following: 

•  Objective  Function  (fun) 

•  Starting  Point  (xo) 

•  Design  Variable  (DV) 

•  Equality  Constraints  ( Aeq.beq ) 

•  Inequality  Constraints  (A,b  and  lb,ub) 

1.  Objective  Function  (OF) 

In  order  to  prove  the  certainty  and  performance  of  the  optimization  program,  in 
the  mass-  error  detection  analysis,  three  types  of  OF’s  were  implemented: 


a.  Eigenvector  Objective  Function 


Taken  from  the  Modal  Assurance  Criterion  concept,  which  quantitatively 
checks  the  correlation  between  the  measured  and  the  calculated  mode  shapes,  its 
arguments  lie  between  0  and  1;  those  elements  with  a  value  close  to  1  indicate  that  the 
two  modes  under  comparison  are  quite  close.  Thus  the  MAC  was  subtracted  from  one  to 
produce  a  result  that  would  indicate  a  better  performance  when  minimized. 
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b.  Eigenvalue  Objective  Function 


Natural  Frequencies  Ratio 


A<p)  = 


i= i 


(7.8) 
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c. 


Eigenvalue  plus  Eigenvector  Objective  Function  Combination 


Eigenvalues  and  eigenvectors  were  selected  as  the  main  core  of  the  OF’s 
since  they  represent  the  modal  response  to  the  natural  properties  of  the  system. 
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2.  Starting  Point  (xo) 

Since  it  was  always  desired  to  iterate  over  the  complete  set  of  elements,  the 
starting  point  consisted  of  a  [42x1]  unit  value  vector. 

3.  Design  Variables  (DV) 

Being  the  eigenvalues  and  eigenvectors  part  of  the  natural  properties  of  the 
system  and  since  these  are  determined  by  the  mass  and  stiffness  of  the  structure,  it  was 
decided  to  use  the  density  (<p)  as  the  design  variable  (DV).  It  was  not  intended  to  modify 
the  elasticity  modulus  (E)  of  the  beam,  such  that  it  could  be  used  in  future  experiments. 
The  design  variable  dv_rho  was  used  under  the  modelo.m  MATLAB  program  to 
assemble  the  finite  element  model  with  the  corresponding  density  perturbation 
introduced. 

4.  Equality  Constraints 

No  equality  constraints  were  used  for  this  experiment,  hence  its  corresponding 
brackets  in  the  matlab  code  were  left  blank. 

5.  Inequality  Constraints 

A  beam  weight  inequality  constraint  was  set  up  in  order  to  bound  the  weight  of 
each  element  by  a  10%,  as  well  as  the  total  amount  of  mass  that  could  be  added  to  the 
entire  beam  by  a  10%. 
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D.  RESULTS  AND  DISCUSSION 


The  three  Objective  Functions  were  minimized  subject  to  the  following 
constraints: 

•  System  totally  unconstrained. 

•  System  with  limited  amount  of  mass  added  to  any  one  element  (10 
percent). 

•  System  with  limited  total  amount  of  mass  added  to  the  entire  beam  (10 
percent). 


The  blue  bars  of  the  following  figures  represent  the  normalized  damage  detection 
distribution  along  the  cantilever  beam  model.  The  top  graphs  of  each  figure  represent  the 
analytical  damage  detection,  while  the  lower  graphs  represent  the  experimental  damage 
detection  distributions  from  the  laboratory.  The  stem  plots  with  the  yellow  circles  on  the 
top  indicate  the  magnitude  and  localization  of  the  actual  errors. 

1.  Modal  Assurance  Criterion  (MAC)  Objective  Function 

Figures  98,  99  and  100  represent  the  analytical  versus  experimental  comparison  of 
the  MAC  objective  function  when  it  was  used  for  the  system  totally  unconstrained,  for  a 
one  element  mass  constrained  system  and  for  the  total  mass  constrained  system, 
respectively. 

The  eigenvector  objective  function  proved  to  work  very  consistently,  since  the 
maximum  value  in  the  table  was  detected  at  element  35  where  the  actual  error  was 
located,  although  there  was  some  error  influence  at  the  neighboring  elements  (34  and  36). 
The  best  damage  detection/  error  prediction  performance  was  found  at  the  total  mass 
constrained  condition,  as  it  is  shown  in  Figure  100,  where  the  error  is  literally  minimized 
and  bounded  at  element  35,  with  no  big  influence  on  the  rest  of  the  elements  of  the 
cantilever  beam. 
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Figure  99.  Analytical-  Experimental  Element  Mass  MAC  O.F. 
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Figure  100.  Analytical-  Experimental  Total  Mass  MAC  O.F. 


2.  Combined  MAC  and  Eigenvalue  Objective  Function 

Figures  101,  102,  and  103  represent  the  analytical  versus  experimental 
comparison  of  the  MAC-  Eigenvalue  objective  function  when  it  was  used  for  the  system 
totally  unconstrained,  for  a  one  element  mass  constrained  system  and  for  the  total  mass 
constrained  system,  respectively. 

This  combined  objective  function  proved  to  work  almost  as  consistent  as  the 
MAC  objective  function  did,  yet  again,  the  more  constrained  the  problem  became  the 
more  isolated  the  error  was  bounded  as  it  can  be  seen  in  Figure  103.  Theoretically 
speaking  the  more  refined  the  objective  function  becomes,  the  better  results  are  going  to 
be  found,  in  this  case  this  combined  objective  function  did  not  perform  better  than  the 
MAC  objective  function  alone,  the  most  likely  reason  is  due  to  the  fact  that  the 
Eigenvalue  objective  function  (equation  7.8)  is  not  accurately  predicting  the  damage 
therefore  is  not  positively  contributing  to  the  combined  objective  function  so  it  can  build 
up  a  better  performance. 
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Figure  101.  Analytical-  Experimental  Unconstrained  M AC/Eigenvalue  O.F. 


Figure  102.  Analytical-  Experimental  Element  Mass  MAC/Eigenvalue  O.F. 
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Figure  103.  Analytical-  Experimental  Total  Mass  MAC/Eigenvalue  O.F. 


3.  Eigenvalue  Objective  Function 

Figures  104,  105,  and  106  represent  the  analytical  versus  experimental 
comparison  of  the  Eigenvalue  objective  function  when  it  was  used  for  the  system  totally 
unconstrained,  for  a  one  element  mass  constrained  system  and  for  the  total  mass 
constrained  system,  respectively. 

The  use  of  the  natural  frequencies  as  objective  function  did  not  provide  a  good 
error  estimate  and  seemed  to  work  as  the  worst  objective  function  out  of  the  three 
selected.  The  main  problem  was  the  exceeded  number  of  iterations  with  the  program 
routine;  hence  this  caused  the  program  to  be  stopped,  before  meeting  the  desired 
tolerance  conditions.  The  best  error  prediction  was  contained  within  the  seventy  percent 
of  the  actual  error,  not  a  very  reliable  situation. 
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Figure  104.  Analytical-  Experimental  Unconstrained  Eigenvalue  O.F. 


Figure  105.  Analytical-  Experimental  Element  Mass  Eigenvalue  O.F. 
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Figure  106.  Analytical-  Experimental  Total  Mass  Eigenvalue  O.F. 


The  objective  function  that  proved  to  work  the  best  was  the  MAC;  the  combined 
MAC-  Eigenvalue  objective  function  did  improve  the  results  from  just  the  natural 
frequencies  objective  function,  but  did  not  get  better  results  with  respect  to  the  MAC 
objective  function  alone. 

The  MAC  objective  function  verified  to  work  more  efficiently  than  the 
eigenvalues  objective  function,  proving  the  simplicity  and  low  cost  of  the  natural 
frequencies  ratio  objective  function,  but  at  the  same  time  it  confirms  the  disadvantages  in 
the  damage  detection  limitation,  because  the  natural  frequencies  can  not  provide  the 
spatial  information  about  structural  damage.  On  the  other  hand,  the  mode  shapes 
objective  function  is  more  sensitive  to  damage,  due  to  the  strain  energy  correlation  at  that 
location  (Jaishi  and  Ren,  2006).  Further  more  the  MAC,  proved  its  usefulness  in  the 
mode  shapes  correlation,  minimizing  the  objective  function  in  question. 
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Theoretically  speaking  the  more  refined  objective  function  should  give  the  best 
damage  detection  results,  but  as  it  was  proved  with  the  combined  MAC-  Eigenvalue 
objective  function,  the  reality  not  always  matches  with  the  theory,  because  it  is  not  only  a 
matter  of  having  a  more  complete/  refined  objective  function,  but  it  is  also  required  to 
assure  that  the  chosen  parameters  assure  reliable  results  by  themselves  in  order  to 
improve  the  results  when  combined  with  other  parameters  and  not  to  weaken  the  actual 
performance  of  the  more  complex/  refined  objective  function.  Even  when  we  had  the  best 
computer  program,  if  the  objective  function  is  not  refined  enough,  the  results  will  not  be 
satisfactorily  found. 
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VIII.  CONCLUSIONS  AND  RECOMMENDATIONS 


Artificial  Boundary  Conditions  have  been  verified  to  work  as  a  method  to  find 
additional  natural  frequencies  by  using  only  one  experimental  data  base,  as  opposed  to 
reconfiguring  the  complete  experiment  such  that  more  than  the  baseline  frequencies 
could  be  found.  Furthermore,  the  main  purpose  of  this  thesis  was  to  accomplish  the 
damage  detection/  error  prediction  on  sensitivity  based  finite  element  models  by  the 
judicious  use  of  the  Artificial  Boundary  Conditions. 

A.  CONCLUSIONS 

1.  The  use  of  Artificial  Boundary  Conditions  via  the  reduced  order  model 
proved  to  be  an  effective  method  to  solve  for  the  spatially  incomplete  real 
world  structures.  Hence  reducing  the  disadvantages  of  the 
underdetermined  problems. 

2.  The  sensitivity  based  updating  governing  equation  proved  to  be  a  very 
useful  tool  in  localizing  the  error/  damage,  when  the  Artificial  Boundary 
Conditions  frequencies  were  added  to  the  baseline  frequencies. 

3.  The  stiffness  sensitivity  distribution  was  directly  correlated  to  the  strain 
energy  distribution  of  the  system;  consequently  the  stiffness  sensitivity  is  a 
function  of  the  flexural  properties  of  the  system  structure. 

4.  The  mass  sensitivity  distribution  was  directly  correlated  to  the  kinetic 
energy  distribution  of  the  system;  as  a  result  the  mass  sensitivity  is  a 
function  of  the  inertia  of  the  system  when  an  excitation  frequency  is 
imposed  over  the  system  structure. 

5.  Damage/  errors  were  very  likely  to  be  detected  at  areas  of  high  sensitivity 
distribution  of  the  structure. 
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6.  The  damage  detection/  error  prediction  was  consistently  improved  when 
the  ABC’s  were  selected  in  a  manner  that  increased  sensitivity  in  the 
desired  region. 

7.  Stiffness  errors  proved  a  very  consistent  behavior.  The  highest  damage 
detection  probability  was  when  the  stiffness  was  the  closest  to  the  fixed 
end  of  the  experimental  cantilever  beam.  The  use  of  the  pinned  ABC  at  a 
DOF  adjacent  to  the  element  error,  consistently  improved  the  damage 
detection  possibilities.  When  the  pinned  ABC  was  applied  to  a  DOF 
further  away  with  respect  to  the  location  of  the  damaged  element,  the  error 
was  less  likely  to  be  detected.  The  use  of  multiple  ABC’s  considerably 
improved  the  error  detection  capability  of  the  FEM  problem. 

8.  Mass  errors  did  not  provide  consistent  evidence  of  damage  detection.  The 
only  conclusion  found  was  that  since  the  mass  errors  are  related  to  the 
kinetic  energy  of  the  system,  the  pinned  ABC  needs  to  be  set  up  such  that 
it  is  able  to  produce  a  considerable  displacement  at  the  location  of  the 
damaged  element.  Pinned  ABC’s  applied  to  mass  based  sensitivity 
systems  tends  to  diverge  the  sensitivity,  hence  any  ABC  applied  at  a  DOF 
adjacent  to  the  mass  damaged  element,  is  less  likely  to  find  the  error 
perturbation.  The  use  of  multiple  ABC’s  did  not  show  an  improvement  as 
far  as  the  error  detection  capability  is  concerned. 

9.  The  optimization  program  considerably  improved  the  mass  error 
detection  performance  of  the  cantilever  beam  model.  The  more  refined  the 
Objective  Function,  the  more  accurate/  minimized  results  are  produced. 

10.  The  modal  to  spatial  transformation  in  the  Optimization  routine  was 
accomplished  by  the  use  of  the  “Simple  Method”,  although  this  method 
proved  to  work,  it  is  believed  that  some  more  accurate  methods  would 
give  a  more  accurate  data  transfer/conversion,  as  far  as  the  softwares 
interface  between  the  MEscopeVES  -  MATLAB  is  concerned. 
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B.  RECOMMENDATIONS 


1.  Find  a  method  of  approach  to  reliably  predict  the  damaged  mass  elements 
over  the  beam  model. 

2.  Find  or  develop  some  interface  routine  to  successfully  pass  directly  the 
data  generated  between  the  software  used  for  the  finite  element  model 
simulation,  such  MEscopeVES,  and  the  software  used  for  the  optimization 
computer  program  (MATLAB). 

3.  Try  more  than  only  one  physical  parameter  of  the  structure  to  be  used  as 
design  variables  in  the  formulation  of  the  Objective  Function,  to  allow 
isolation  of  all  modeling  errors  while  not  requiring  excessive  iterations, 
hence  obtaining  an  improved  optimization  result. 
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APPENDIX 


A.  ABCRUNTHRURLA 


ABCr  u nTHRU  r I  a .  m 


%  This  program  calculates  the  condition  number  of  the  following 
%  sensitivity  matrices  used  to  calculate  the  DV  (error  prediction). 
%  1)  Base  system  only  5  modes  (underdeter  mi  ned) 

%  2)  ABC  system  10  modes 

%  3 )  Base  system  5  modes  +  5  modes  from  ABC  system 
%  The  last  system  is  calculated  3  times.  Once  for  modes  1-5, 

%  another  for  modes  6-10,  and  again  for  modes  11-15. 

% 

%  This  program  is  called  from  Build2Beams.m. 

% 

%  Written  by  Constance  R  S  Fernandez,  Spring  2004 
%  Updated  by  John  R  Mentzer,  Spring  2007 


%l  NPUTS 


%  i  c  n  t  o  s  e  t 
%  T  sens  tot 
%  v  e  c  t  _  I  a  m_  t  o  t 

%  OUTPUTS 


%  a  a 
%  dv  a 
%  d v" b 
%  dv"c 
%  i  n  t  e  r  v  e  I 
%  mode 
%  s  t  a  r  t  mo  d  e 
%  ten 

%  dv  cal  ABC  -  mat  r  i  x 
%  d  v “ c  a  I  ABCt  en  -  mat  r  i  x 
%dv"cal  BasePlus  -  matrix 
%  c o n d  ABC  -  mat  r  i  x 
%cond"ABCten  -  matrix 
%  cond“basePI  us  -  ma t  r  i  x 

% - I  Nl  Tl  ALI  ZATI  ON . 


a  be  =0; 
ten  =  1; 
i  n  t  e  r  v  e  I  =  1; 
i  nt  abc  =  1; 

f  o  r  “  a  a  =  1 :  i  c  n  t  oset  +1  %  number  of  conditions  (base  +  ABC) 
a_c  =  1;  %  reinitialize  for  each  ABC  system 

for  mode  =  1:3  %  3  sets  of  modes  per  ABC  system  (10  element  beam) 
%  (  mo  d  e  s  1:5,  6:10,  11:15) 
start  mode  =  abc  +  a_c; 

dv  a  =  [start  mode:  st  art  mode+4] ;  %  modes 

dv"dl  =  [ t  en:  t  en] ;  %  mode  1 

d v ” d 2  =  [  t  en+1:  t  en+1] ;  %  mode  2 

d v “ d 3  =  [  t  en+2:  t  en+2] ;  %  mode  3 

d v ” d 4  =  [  t  en+3:  t  en+3] ;  %  mode  4 

d v ” d 5  =  [  t  en+4:  t  en+4] ;  %  mode  5 


dv  c  =  [ten: t  e  n  +9 ] ;  %  the  first  10  modes  of  each  ABC  system 
%  [  mo  d  e  s  1  - 1 0  o  n  I  y ) 

if  dv  a  ==  [  1:  .  25*si  ze(T  sens  tot, 2)];  %  if  sensitivity  matrix 
%fias  only  5  rows  than  the" modes  used  are  all  5,  else  modes 
%used  are  first  five  (base)  and  a  set  of  selected  5  modes 
%of  ABC  system  for  a  total  of  10  modes, 
dv  b  =  [ dv  a ] ; 
else" 
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dv  b  =  [ 1: .  25*si  ze(T  sens  t  o  t , 2 ) ,  dv  a]; 
end 

% —  Base  System  onl  y  —  (underdetermi  ned) 

%  save  DV  calculates  of  as  matrix  for  plotting 

dv  cal  A  B  C ( :  ,  i  ntervel  )  =T  sens  tot(dv  a ,  :  )  \  v  e  c  t  lamtotfdv  a);  hold  on 
%  condition  number  of  Sensitivity  matrix  used  in'DV  cal. 
cond  A  B  C  ( i  ntervel  ,  1)  =  cond(T  sens  tot(dv  a,:)); 

%  cond  A B C ( i  nt  er  vel ,  1)  = 

c  o  n  d ( ( T  sens  tot(dv  a ,:))*(  T  sens  tot(dv  a,  :))');  %(  T  t  r  a  n  s )  (  T ) 


% 

%mo  d  e  1 

% 

dv  cal 

ABC1  =  T_ 

s  ens_ 

t  o  t  ( d  v  _ 

dl, : ) \ vect 

_  1  a  m_ 

t  o  t  ( d  v  _ 

dl) 

% 

%mo  d  e  2  " 

% 

dv  cal 

ABC2  =  T_ 

s  e  n  s  _ 

t  o  t  ( d  v  _ 

d2, : ) \ vect 

_  1  a  m_ 

t  o  t  f  d  v  _ 

d  2 ) 

% 

%mode  3" 

% 

d  v  cal 

A  B  C  3  =  T_ 

sens_ 

t  o  t  ( d  v  _ 

d3, : ) \ vect 

_  1  a  m_ 

t  o  t  ( d  v  _ 

d  3 ) 

% 

%mode  4" 

% 

dv  cal 

A  B  C  4  =  T_ 

s  ens_ 

t  o  t  ( d  v  _ 

d  4 , : ) \ vect 

_  1  a  m_ 

t  o  t  ( d  v  _ 

d  4 } 

% 

%mo  d  e  5  " 

% 

dv  cal 

A  B  C  5  =  T_ 

s  ens_ 

t  o  t  ( d  v  _ 

d5, : ) \ vect 

_  1  a  m_ 

t  o  t  f  d  v  _ 

d  5 ) 

%  ---Base  +  5  modes  of  ABC  system  - 

dv  cal  BasePI  us(  :  ,  i  ntervel  )  =T  sens  tot(dv  b,:)\vect  lamtotfdv  b); 
cond  basePI  usfi  ntervel  ,  1)  =  condfT  sens  tot[dv  b,:));~ 

%  cond  basePI  usfi  ntervel,  1)  = 

condff (T_sens“tot(dv_b, : ))*(T_sens_tot(dv_b,  :  ))  1  ));  %(  Tt  r  a  ns )  (  T) 

i  ntervel  =  i  n  t  e  r  v  e  I  +1; 
a _ c  =  a _ c  +  5 ;  %  f  i  v e  mo d e s  used  at  a  t  i  me 

end  %  11  mo  d  e 11  loop 

%  ---ABC  system  only  - 


dv  cal  ABCtenf  :  ,  i  nt  abc)  =T  sens  totfdv  c,:)\vect  lamtotfdv  c); 
cond  ABCtenf i  n t  a b c 7 1 )  =  condfT  sens  tot[dv  c,:));~ 

%cond  ABCtenfint  abc,l)  = 

c  ondf ( ( T_s  ens~t  ot ( dv_c , : I ) * ( T_s  ens_t  ot ( dv_  c ,  :  ) )  1  ) ) ;  %(  Tt  r  a  ns )  (  T) 
i  nt  abc  =  i  nt  abc  +  1; 

abc"=  abc  +  39;  %  advances  to  next  ABC  system,  must  change  number  "19" 
%to  reflect  the  number  of  DOF  in  beam.  This  beam  had  10  elements  thus 
%  19  DOF. 

ten  =  ten  +  39;  %  advances  to  next  ABC  system 
end  %  11  a  a 11  loop 


END  ABCrunTHRU  rl  a.  m 


118 


B. 


ADDLUMPMASS  RLA 


%  ********************  AddLumpMass  rla.m  *********************** 

%  This  script  constructs  a  vector  of  I  u mp e d  ma s s e s 
%  which  is  added  to  the  diagonal  of  the  BeamX  mass  matrix. 

%  Mass  added  to  [  mx  ]  in  A  s  s  e  mb  I  e  2  B  e  a  ms .  m 

% 

%  Wr  i  1 1  en  by  Prof  J  .  H.  Gor  di  s 

%  I  n  p  u  t  s  needed 

% . 

%  n  u  m  e  I  e  me  n  t  s 
%  mx 


%  Out  put  s 

% . 

%  ma  s  s  d  i  a  g 
%  mass" node 
%  ma  s  s "  c  o  o  r  d 
%  mass" DOF 
%  mx  -"updated 


di  sp( 
di  sp( 
di  sp( 
di  sp( 
di  sp( 


di  sp( 


****  Lumpedmassadditiontobeams  ****') 

********  S?1**1****  *  *  * *************************  *****?*  *******  j 


%  i  n  i  t  a  I  i  z  e 

if  exist(' mass_diag')  ==  0;  %  define  and  apply  lumped  mass  vector. 


add  mass  =  1  n1  ; 

a d d “ ma s s  =  input!1  Add  lumped  masses  to  BeamX  ?  (y/n)  Vs1); 


%  Initialize  vector  to  add  to  [mx]  diagonal. 


ma s s _ d i  a g  =  zeros!  2*(  num_el  ements+1) ,  1) ; 


whi  I  e  a d d _ ma s s  ==  1  y1  ; 

mass_node  =  input!1  Node  number  for  lumped  mass  ?  1  ) ; 


mass_coord  =  input!1  Translation  or  Rotation  for  lumped  mass  ?  (t/r)  1 


if  mass  coord  ==  't1;  %  Tr  a  ns  I  at  i  ona  I  lumped  mass 
mass  DOF  =  2  *  mass  node  - 
el  seif  mass  coord  ==  'rr;  %  rotational  lumped  mass 
mass  DOF  =2  *  mass  node; 


mass  d  i  a  g  (  mass  DOF)  =  input!1  Enter  value  of  mass/inertia  (in  "lbf-secA2/ 
% "  p  u  t  s  lumped  mass  on  correct  DOF 
add  mass  =  input!1  Add  another  lumped  mass  ?  (y/n)  Vs1); 

%  can  continue  adding  mass  until  'n'  is  entered 

end;  %  End  whi  I  e  loop 

end;  %  End  e  x  i  s  t  ( 1  ma  s  s  _  d  i  a  g 1  ) 

mx  =  mx  +  d i  a g (  ma s s _ d i  a g ) ;  %  A d d  lumped  ma s s e s  to  [  mx ]  : 

0^  ********************  END  ADDLUMPMASS  rla.m  *********************** 
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c. 


ASSEMBLESENS  RLA 


%  ********************  A  s  s  e  mb  I  e  S  e  n  s  r  I  a .  m  *********************** 

% 

%  This  program  asse mb les  the  total  sensitivity  matrix,  T  sens  tot  and 
%  total  lam  vector,  vect  lam  tot  and  asse mb  les  the  relative  frequency 
%  error  between  the  natural  frequencies  of  Beam  A  and  BeamX 
%  Written  by  Constance  R  S  Fernandez,  Spring  2004 
%  Updated  by  John  R  Mentzer,  Spring  2007 


%  I  NP UTS 

% . 

%  v  e  c  t  I  a  mx  o  s  e  t 
%  v  e  c  t "  I  a  m  oset 
%  v  e  c  t "  I  a  m" 

%  T  sens  oset 
%  T"sens" 

%  " 

%  OUTPUTS 

% . 

%  vect  OSET 
%  v  e  c  t  I  a  m  t  o  t 
%  T  sens  tot 
%  f  r  e q  OSET 
%  f  r  e q" OS  ETx 
%  r el  _f r  e q ERROR 

vect  OSET  =  vect  lamx  oset  -  vect  lam  oset; 

%  I  a  mx  from  actual  beam  with  error  oset,  lam  from  FE  model  oset 
%  Creating  a  vector  of  lam  differences  calculated  ( Lx-  La) 
if  vect  OSET  ==  0; 

%when  vector  is  empty  at  first,  the  total  vector  is  equal  to  the 
%  lam  vector  of  Beam  A,  i .  e . ,  the  first  19  values  of  vect  lam  tot  are 
%  t  h  e  natural  freq  squared  (radA2/secA2)  of  BeamA 

vect  lam  tot  =  vect  lam; 

else 

vect  lam  tot  =  c  a  t  ( 1 ,  vect  lam,  vect  OSET); 
end 


if  T_s  e  n  s  _  o  s  et  ==  0; 

T  sens  tot  =  T  sens; 

else 

T  sens  tot  =  c  a  t  ( 1 ,  T  sens,  T  sens  oset); 
end 

freq  OSET  =  s q r t ( a bs ( vec t  OSET))/2/pi; 

%  Natural  frequency  vector  of  BeamA  in  Hz 
freq  OSETx  =  s  q  r  t  ( a  bs  (  vec  t  lamx  oset))/2/pi; 

%  Natural  frequency  vect  or  "of  Beam  X  in  Hz 
r  el  freqERROR  =  freq  OSET. /freq  OSETx*100; 

%  Relative  error  between  Beam  A  OSET  natural  freq  and  BeamX  OSET  natural  Freq. 

o/0  ********************  end  Asse  mb  I  eSens  rla.m  *********************** 
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D. 


BEAMPROPERTIES  RLA 


%  ********************  BeamProperti  es  rla.m  *********************** 

%  This  is  the  "props  file"  to  load  nominal  beam  data. 

%  This  program  is  cafled  by  BeamA  Prompt  crs  to  provide  beam  properties 
%in  order  to  build  BeamA. 

%  Thi  s  program  was  written  by  Constance  Fernandez,  Spring  2004 
%  This  progam  mo dified  by  John  R  Mentzer,  Spring  2007 
%  Out  put  s 

% . 

%  depth 
%  wi  dt  h 
%  E 
%  r  h  o 

%  t  o  t  a  I  length 
%  n  u  m  e  F  e  me  n  t  s 
%  n o mF  n a  I  El 
%  n  o  mi  n  a  I  ~  a  r  e  a 
%  n  o  mi  n  a  I "  d  e  n  s  i  t  y 

%  Following  are  actual  measurements  from  experimental  set-up  cantilever 
%  beam. 

depth  =  .504;%  in  z-dir  (inches) 
width  =  1.506;  %in  y-dir  (inches) 

E  =  1 0  e  6 ; 

%E  =  1.  6 5 e 6 ;  %  I  bf  /  sec A2- i  n  ( 1 0 e 6 - ksi  ) 

%(  1  bf  /  i  n  A2  =  6  8  9  4.  7  6  P a )  -  >  E (  I  bf  /  i  n  A2 )  =  ()  Pa/  6  8  9  4.  7  6 

r ho  =0. 1  1  0  4  6  0  9  3  4;  %0.  0  9  8;  %  I  bf / i  nA3 

%T6  temper  alloys  require  a  3  5  -  ksi  tensile  strength,  3  0  -  ksi  yield 
%  strength  and  a  1 0  e  6  -  ksi  elastic  modulus.  Alloy  6061-T6  has  1.0 

%  pet  magnesium,  0.6  pet  silicon,  0.3  pet  copper  and  0.2  pet  chromium. 

%  It  has  a  4  5  -  ksi  tensile  strength  and  3  5  -  ksi  yield  strength. 1  The 
%  machinability  of  aluminum  alloys  are  high  (300)  compared  to  titanium 
%  (40).  Aluminum  alloys  can  easily  be  bent  and  provide  easy  loading  and 
%  unloading  of  parts.  Also,  aluminum  is  a  highly  conductive  metal 
%  c  o  mp  a  r  e  d  t  o  t  i  t  a  n  i  u  m. 


%  a  I  I  measurement  of  distance  are  in  inches 


total  length 
num  efements 
n o mf  n a  I  El 
n  o  mi  n  a  I"  a  r  e  a 


=  42 


=  20; 

=  ( wi  dt  h  *  dept  h A3  /  12)  * 
=  depth  *  wi  dt  h;  %  i  n  A2 


nomi  nal  _densi  ty  =  r  h  o ;  %  I  b  f  /  i  n  A3 


E; 


^  ********************  END  B  e  a  mP  ropert  i  es  rla.m  *********************** 
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E. 


BEAMSENSITIVITY  RLA 


%  ********************  BeamSensitivity  rla.m  ************** 
%  Revision  history: 

% 

%  Ver.  1.0:  4/4/95  Basic  frequency  sensitivities 
%  Updated:  Spring  2004  Constance  R  S  Fernandez 


%  Program  Des  c  r  i  pt  i  on : 

% . 

% 

%  This  program  calculates  mode  frequency  sensitivities  as  given  by  the 
%  equation, 

% 

%  1wA2  HI  k]  1  [  m] 

%  --  =  {P}1  *  [  . wA2  ]  *  {P} 

%  1DV  H  DV  H  DV 

% 

%  where:  w  =  natural  frequency 
%  {P}  =  associated  mode  shape 

%  DV  =  desi  gn  var i  abl  e 

% 

%  The  right  side  of  the  above  equation  is  considered  the  addition  of 
%  the  sensitivity  matrix  wrt  El  and  the  sensitivity  matrix  wrt  mass. 

% 

%  The  program  calculates  the  stiffness  and  mass  matrix  partials  by  finite 
%  differences.  That  is,  for  example,  the  [k]  matrix  is  assembled  twice, 

%  once  in  for  the  nominal  beam  parameters,  and  a  second  matrix  is 
%  assembled  based  on  a  small  perturbation  of  element  mass  and/or  El. 

% 

%  This  program  makes  use  of  the  beam  data  created  by  the  program 
%  "  B  u  i  I  d  2  B  e  a  ms .  m.  11 

%  1)  Resets  BeamX  mass  and  El  data  to  be  the  same  as  BeamA  data. 

%  2)  Enters  a  loop  to  create  a  sensitivity  matrix  (T) 

%  I  n  I  o  o  p 

%  a)  A  small  perturbation  of  mass  is  added  to  BeamX  on  ele.l. 

%  b )  The  mass  ma  t  r  i  x  is  a  s  s  e  mb  I  e  d  for  this  ma  s  s  ■  p  e  r  t  u  r  b  e  d  beam, 

%  and  the  mass  matrix  partial  derivative  is  calculated  as 

% 

%  H[m]  [  m  p  e  r  t  u  r  b  ]  ■  [  m  b  a  s  e  I  i  n  e  ] 

%  =  - 
%  H DV  H DV 

% 

%  c)  The  first  column  of  Sensitivity  matrix  is  calculated  using: 

% 

%  sens  mass  =  [phi  base]'*(-lam  base)*m  delta*[phi  base] 

% 

%  (Note:  A  column  of  T  corresponds  to  the  respective  element  on  beam.) 

% 

%  d)  T  loop  starts  again  but  with  the  small  change  on  element  2. 

%  e)  Difference  calculated  and  second  column  of  T  is  calculated. 

% 

%  3)  loop  continues  until  all  columns  of  T  are  calculated,  corresponding 

%  to  a  small  change  added  to  the  respective  element. 

% 

%  4)  The  procedure  is  identical  for  stiffness  sensitivity. 

%  Except: 

%  s  e  n  s  _  E I  =  [phi  base]'*kdelta*[phi  base] 

% 

%  5)  Combine  the  two  T1  s  into  one  complete  sensitivity  matrix  : 

%  T  sens  =  [  s  e  n  s  ma  s  s ,  s  e  n  s  El]. 

%  In  words,  the  first  "set  of  cofumns  (equal  to  number  of  elements)  of 

%  the  combined  matrix  is  the  sensitivity  mass  wrt  mass  changes  and 

%  the  last  half  is  wrt  El  changes. 

% 


0/0  **************************************************************************** 

%  I  n  p  u  t  s  Needed: 

% . 
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%  ma  s  s  Ibis 
%  E I  Ibis 
%  e I  e me n t  mass 
%  e  I  e  me  n  t  "  E I 
%  n  u  m  r  b  m" 

%  n  u  m"  e  I  e  me  n  t  s 

%  I  a  mx  ( e x  p e  r  i  me  n  1 1  y  me  a  s  u  r  e  d  n  a  t  f  r  e  q  of  beam) 

%  P  r  o  g  r  a  ms  called 

% . 

%  A  s  s  e  mb  I  e  2  B  e  a  ms  c  r  s ; 

%  BoundaryCondi tfons  crs; 

%  f  mo  d  e  s 

%  fOset_from_Aset 

%  Out  put  s : 

% . 

%  num  modes 
%  I  a  m"  b  a  s  e 
%  p  h  i "  b  a  s  e 
%  sens  mass 
%  sens" El 
%  T  sens 
%  dv  cal 
%  f  r  e  q  base 
%  v e c t ” I  am 

%  Start  code: 

%  . 

f  o  r  ma  t  long; 


I  Nl  T I  ALI  ZAT I  ON 


mass  change  =  1;  %  Percent  mass  change  for  sensitivity  calculation. 
El_cbange  =  1;  %  Percent  El  change  for  sensitivity  calculation. 

element  mass  orig  =  element  mass;  %  Copy  properties  to  retain  them, 
el  ement”EI  _ori  g  =  elemental; 

%  Prompt  for  number  of  mode  frequencies  to  generate  sensitivities  for: 

% . 


%num  modes  =  input)1  Enter  number  of  modes  for  sensitivity  calculations  >>  ') 
num  modes  =39;  %  3  9  is  maximumnumber  of  modes  in  a  20  element  cantilever 
%beam.  This  is  hard  coded  for  quicker  run  of  simulation. 

start_mode  =  num_rbm+  1;  %  Skip  the  rigid  body  modes. 

%**************  MASS  SENSITIVITY  CALCULATION  LOOP  ******************* 

sens_mass  =  0;  %  initialize  Mass  based  sensitivity  matrix 

if  ma  s  s  _  I  b  I  s  ~=  0;  %  From  Beam_XPrompt  as  user  inputs 

for  i  c  n  t  _  d  v  =  1:  num_el  ements;  %  loop  to  create  sensitivity  matrix 

%  Resetting  BeamX  properties  to  BeamA  properties 
el  ement  mass) :  ,  2)  =  e I  e me n t _ ma s s ( :  ,  1 ) ; 

%  each  element,  one  at  a  time  will  have  a  change  in  mass 
el  ementmass)  i  cnt  dv,2)  =  element  mass)  i  cntdv,  2)  *  ... 

(  1  +  mass_change/ 100) ; 

Assembl  e2Beams_crs;  %  Run  script  to  assemble  beams. 

BoundaryCondi  ti  ons_crs;  %  Apply  boundary  conditions. 

[la  m_  base,  phi  _  b  a  s  e  ]  =fModes(ka,ma); 

%  ma  is  the  basebeam  wi  t  hout  changes, 

%  mx  is  beam  with  small  change  in  mass  added  (  mass_c  hange) 

%  Form  ma  s  s  derivative  ma  t  r  i  c  e  s : 
m_  d  e  I  t  a  =  (  mx  -  ma )  /  (  ma  s  s  _  c  h  a  n  g  e  / 1 0  0 ) ; 


%  converts  precent  change  into  decimal  a  mount 


%  Mode  f  r  e q  sens  loop: 
end  mode  =  start  mode  +  (num  modes  ■  1); 
r  o  w"  n  u  m  =  0 ;  %  initialize  loop 
for"icnt  modes  =  start  mode:end  mode; 
row  num  =  row  num  +1; 

sens  mass  ( row- num,  icnt  dv)  =  phi  base(:,icnt  modes)1  *... 
[-Iambase(icnt  modes)  *  m  delta  j  *... 
phi  b  a  s  e  ( :  ,  icnt"  modes); 

%  definTtion  of  mass  sensitivity  ma  t  r  i  x 
end;  %  f  o  r  11  i  c  n  t  _  mo  d  e  s 11  inner  loop 

end;  %  End  "for  i  c  n  t  _  d  v 11  outer  loop  for  sensitivity  calculations 
end;  %  End  11  i  f  ma s s _ I  b I  s  ~=  0" 


0^  El  SENSI T I  VI TY  CAL CUL ATI  ON  LOOP 

sens  El  =0;  %  initialize  El  based  sensitivity  matrix 

if  El  _  I  bis  ~=  0;  %  From  Beam_XPrompt  as  user  inputs 

for  i  c  n  t  _  d  v  =  1:  num_el  ements;  %  loop  to  create  sensitivity  matrix 

%  Resetting  BeamX  properties  to  BeamA  properties 
element_EI(:,2)  =  element_EI(:,l); 

%  each  element,  one  at  a  time  will  have  a  change  in  El 
element  Elficnt  d  v ,  2 )  =  element  Elficnt  d  v ,  2 )  *  ... 

( 1  +  (El  _  c  h  a  n  g  e/ 100)  ) ;  “ 

Assembl  e2Beams_r  I  a;  %  Run  script  to  assemble  beams. 

BoundaryCondi  ti  ons_rl  a;  %  Apply  boundary  conditions. 

[la  m_  base,  phi  _  b  a  s  e  ]  =fModes(ka,ma); 


%ka  is  the  basebeam  without  changes, 

%kx  is  beam  with  small  sensitivity  added 

%  F  o  r  m  E I  derivative  ma  t  r  i  c  e  s : 
k  delta  =  (kx  ■  k  a )  /  ( El  change/100); 

%" converts  precent  change  to  decimal  value 

%  Mode  f  r  e q  sens  loop: 

end  mode  =  start  mode  +  (num  modes  ■  1); 

r  o  w"  n  u  m  =  0 ; 

for  icnt  modes  =  start  mo  d  e :  e  n  d  mo  d  e  ; 
row  hum  =  row  num  +1; 

sens  El  ( row  num,  icnt  dv)  =  phi  b  a  s  e  ( :  ,  i  c  n  t  modes)1  *... 

k  delta-*  phi  b  a  s  e ( : , i  c  n  t  modes); 

%d  e  f  i  n  i  t  i  on  of  El  "sensi  ti  vi  ty"matri  x 
end;  %  f  o  r  11  i  c  n  t  _  mo  d  e  s 11  inner  loop 

end;  %  End  "for  icnt_dv"  outer  loop  for  sensitivity  calculations 

end;  %  if  EM  bl  s  =0; 


%  Copy  element  El  and  mass  properties  back  into  arrays: 
element  El  =  e  I  e  me  n  t  El  orig; 
e  I  e  me  n  t "  ma  s  s  =  element"mass_orig; 

%  cleans  up  workspace  by  clears  all  unimportant  parameters 
clear  element  El  _  o  r  i  g  element  mass  orig  end  mode  start  mode 
clear  row  num- icnt  modes  El  change-k  delta  mass  change"m  delta 
clear  ka  l<x  ma  mx  fcnt_dv 

%  A  s  s  e  mb  I  e  s  complete  sensitivity  ma  t  r  i  x 
if  sens_mass  ==  0&  s e n s _ E I  ~=0 ; 

T  sens  =  sens  El;  %  resultant  sensitivity  matrix  is  equal 
%"to  El  sensitivity  matrix  with  only  El  changes  if  no  inputs 
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%  are  given  for  mass  changes 

elseif  sens  mass  ~=  0  &  sens  El  ==  0; 

T  sens  =  sens  mass;  %  resultant  sensitivity  matrix  is  equal 
%"to  Mass  sensitivity  matrix  with  only  mass  changes  if  no  inputs 
%  are  given  for  El  changes 


else 

T  sens  =  c  a  t  ( 2 ,  sens  mass, sens  E I  ) ;  %  else  the  resultant  sensitivity 
%ma trix  is  teh  combination  of  mass  sensitivity  matrix  in  first  set 
%of  columns  and  El  sensitivity  matrix  in  the  last  set  of  columns 


end 

freqx  =  sqrt(lamx)/2/pi; 

%freq  base  =  sqrtjlam  base)/2/pi; 

%  NOTE:  %  lamx  =  experiment  measured  natural  freq  of  beam 
v  e  c  t  _  I  a  m  =  ( I  a  mx  ( 1 :  n  u  m_  mo  d  e  s )  - 1  a  m_  b  a  s  e  ( 1 :  n  u  m_  mo  d  e  s ) ) ; 

o/0  ********************  end  BeamSensitivity  rla.m  *********************** 
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F. 


BE AMSEN SITIVIT YOSET  RLA 


BeamSensi  t  i  vi  t  y OS E T _ r  I  a.  m 


%  Revision  history: 

% . 


%  Ver.  1.0:  4/4/95  Basic  frequency  sensitivities 
%  Updated:  Spring  2004  Constance  R  S  Fernandez  to  create  resulting 
%  sensitivity  matrix  using  lam  vector  of  multiple  ABC  systems 


%  Program  Des  c  r  i  pt  i  o n : 

% . 

% 

%  This  program  calculates  mode  frequency  sensitivities  as  given  by  the 
%  equation, 

% 

%  HwA2  1  [  k ]  1  [  m] 

%  ■■  =  { P } 1  *  [  . wA2  -  -  -  ]  *  { P } 

%  1DV  HDV  ID V 

% 

%  where:  w  =  natural  frequency 
%  {P}  =  associated  mode  shape 

%  DV  =  design  variable 

% 

%  The  right  side  of  the  above  equation  is  considered  the  addition  of 
%  the  sensitivity  matrix  wrt  mass  and/or  El. 

% 

%  The  program  calculates  the  stiffness  and  mass  matrix  partials  by  finite 
%  differences.  That  is,  for  example,  the  [k]  matrix  is  assembled  twice, 

%  once  in  for  the  nominal  beam  parameters,  and  a  second  matrix  is 
%  assembled  based  on  a  small  perturbation  of  element  mass  and/or  El. 

% 

%  This  program  makes  use  of  the  beam  data  created  by  the  program 
%  "Bui  I  d  2  B  e  a  ms .  m.  11 

%  1)  Resets  BeamX  mass  and  El  data  to  be  the  same  as  BeamA  data. 

%  2)  Enters  a  loop  to  create  a  sensitivity  matrix  (T) 

%  I  n  I  o  o  p 

%  a)  A  small  perturbation  of  mass  is  added  to  BeamX  on  ele.l. 

%  b )  The  mass  ma  t  r  i  x  is  a  s  s  e  mb  I  e  d  for  this  ma  s  s  -  p  e  r  t  u  r  b  e  d  beam, 

%  and  the  matrix  partial  derivative  is  calculated  as 

% 

%  H  [  m]  [  m  p  e  r  t  u  r  b  ]  -  [  m  b  a  s  e  I  i  n  e  ] 

%  =  "  - 
%  HDV  HD V 

% 

%  c)  The  first  column  of  Sensitivity  matrix  is  calculated  using: 

% 

%  sens  mass  =  [phi  base]l*(-lambase)*mdelta*[phi  base] 

%  ~  " 

%  (Note:  A  column  of  T  corresponds  to  the  respective  element  on  beam.) 

% 

%  d)  T  loop  starts  again  but  with  the  small  change  on  element  2. 

%  e)  Difference  calculated  and  second  column  of  T  is  calculated. 

% 

%  3)  loop  continues  until  all  columns  of  T  are  calculated,  corresponding 

%  to  a  small  change  added  to  the  respective  element. 

% 

%  4)  The  procedure  is  identical  for  stiffness  sensitivity. 

%  Except: 

%  sens  El  =  [phi  base] 1  *  k  del  t  a  *  [  phi  base] 

%  _  _  _ 

%  5 )  C  o  mb  i  n  e  the  two  T 1  s  into  one  complete  sensitivity  matrix  : 

%  T  sens  =  [  s  e  n  s  ma  s  s ,  s  e  n  s  El]. 

%  In  other  words,  the"first  set"of  columns  (equal  to  number  of  elements)of 

%  the  combined  matrix  is  the  sensitivity  with  respect  to  (wrt)  mass 
%  changes  and  the  last  half  is  wrt  El  changes. 

% 


%  I  n  p  u  t  s  Needed: 

% . 

%  ma  s  s  Ibis 
%  E I  I  E I  s 
%  e I  e me n t  mass 
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%  el  ement  El 
%  n  u  m  r  b  m" 

%  mx  beam 
%  kx'beam 

%  P  r  o  g  r  a  ms  called 

% . 

%  A  s  s  e  mb  I  e  2  B  e  a  ms  c  r  s ; 

%  BoundaryCondi tfons  crs; 
%  f  mo  d  e  s 

%  f  Oset  from  Aset 
%  di  spl  acment  PI  o  t  _  OS  E  T 

%  Out  put  s : 

% . 

%  num  modes  0 
%  sens  mas s 0,  sens  El  0 
%  num  rbmOSET 
%  T  sens  oset 
%  dFspX  tot,  dispA  tot 
%  a  c  c  e  I  "  p  I  o  t 
%  o  s  e  t  choice 
%  a  c  c  e  F  o  me  t  e  r 
%  BC,  BCose,  BCOSET 
%  r  e  ma  i  n  d  o  f 

%  ASETtot,  OSET,  OSETtot 
%  i  n  c  t  sens 

%  mass"change,  El  change 
%  phi  XPLOT,  phi APLOT 
%  maO  base,  kaObase 
%  mxO,  kxO 
%  p  I  o  t  mx ,  p  I  o  t  k  x 
%  I  amaOSET,  phi  a  OSET 
%  I  amxOSET,  phi  zOSET 
%  I  a  mx  p  I  o  t ,  p  h  i  x  p  I  o  t 
%  f  a  0 

%  m_  d  e  I  t  a  0 ,  k  _  d  e  I  t  a  0 

%  Start  code: 

%  . 


I  Nl  T I  ALI  ZAT I  ON 

T  s ens  os et  =  [ ] ; 
vec  t  lam  oset  =  [ ] ; 
d  i  s  p  X  tot  =  [  ] ; 
di spA"tot  =  [  ] ; 
i  net  sens  =  0; 
icnt"oset  =  0; 

BCOSET  =  z er  os ( ndof ,  ndof ) ; 

OSETtot  =  z  er  os ( ndof ,  ndof ) ; 

ASETtot  =  z  er  os ( ndof ,  ndof); 

oset  choice  =  1  n1  ; 

ndof"%  prints  ndof  to  screen  for  user's  reference 

accelometer  =  [  3 :  2 :  n  d  o  f  ] ;  %This  needs  to  change  if  you  do  not  use  a  clamped  free 
beam! 

%a  c  c  e  I  o  me  t  e  r  =  [3  5  7  9  11  13  15  17  19  2  1  ];  %  use  this  line  for  convenience 

%  in  quicker  calculation  loops  or  use  the  next  2  lines. 

%ac  c  el  o  meter  =  input  ('On  what  nodes  are  the  accel  .  I  ocated'  ,  1  s'  );  %requests 

%  u  s  e  r  t  o  i  n  p  u  t  I  o  c  a  t  i  o  n  s  o  f  a  c  c  e  I  o  me  t  e  r  s 

%a  c  c  e  I  o  me  t  e  r  =  e  v  a  I  ( [ '  [  '  ,  a  c  c  e  I  o  me  t  e  r ,  '  ]  '  ] ) ;  %  converts  string  to  vector  of 

%  I  a  be  I  s 

%  saved  for  plotting  accelometers  in  correct  position, 
accel  plot  =  f  I  oor  ( accel  omet  er  /  2) +1; 

%  graphical  representation  of  accelometer  locations 

oset  choice  =  inputf'  Do  you  want  to  use  ASET  and  OSET  ( y / n ) ?  ' , ' s '  ) ; 

%  user  can  choose  to  use  ASET  and  OSET  calculations 
whi  I  e  oset  choice  ~=  '  n '  ; 

ient  oset  =  ient  oset  + 1 ; 

dispC  I  ocati  ons-of  Accel.')  %  displays  "location  of  Accel"  on  screen 
accelometer  %  displays  the  vector  of  location  of  Accel  for  user's 
%  reference  in  choosing  which  DOF  are  to  be  pinned 
dispC  '  ) 

di  s  p( 1  Enter  pinned  DOF  I  abel  ( s )  1  ) 
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d  i  s  p  (  1  Use  MATLAB  vector  f  o  r  ma  t  >  1  3  5:7  9  ') 

pinned  =  i  nput  ( 1  >>  1  ,  1  s 1  ) ; 

pinned  =  eval  (['  ['  ,  pi  nned,  1  ]'  ]);  %  Converts  string  to  vector  of  labels 
BCficnt  o  s  e  t ,  :  )  =  [  1,  2,  pi  nned];  %  Boundry  conditions  for  cantilever  beam, 

%  change  line  if  different  beam  is  used 

BCoset  =fOset_from_Aset(ndof,  BC(icnt_oset,:));  %  gets  OSET  from  ASET 

B  C  OS  E  T ( i  cntoset ,  1 : I  engt h( BCos et ) )  =  BCoset;  %  copies  OSET  into  another 
%  vector  to"be  used 

%  loop  to  find  r  e  ma  i  n  i  n  g  a  c  c  e  I  o  me  t  e  r  s . 
for  icnt  =  1  :  I  engt  h( pi  nned) ; 

remai  n(i  cnt,  :  )  =  fi  nd(  pi  nned(  i  cnt)  ==  accelometer); 
end 

remaindof  =  f  Os  e  t  from  Aset((l  ength(accel  ometer)),  remain); 

%  r  emai  ni  ng  unpinned  DOF 

%  r  e  ma  i  n  i  n  g  a  c  c  e  I  o  me  t  e  r  s 
aset  =  accelometer(remaindof); 

ASETt ot ( i  c nt _oset ,  1 : I  e n gt h ( a s et ) )  =  aset; 

%  o  mi  1 1  e  d  oset,  i.e.  pinned  accelometers 
OSET  =  f  Os  e  t  f r om  Aset  ( ndof ,  aset); 

OSETtot(icnt"oset7  1:  length!  OSET))  =  OSET;  %  copies  OSET  into  another 
%  vector  to  be  used  i  n====== 

i  nct_sens  =  i  nct_sens+l; 

mass  change  =  1;  %  Percentage  mass  change  for  sensitivity  calculation. 

El_cbange  =  1;  %  Percentage  El  change  for  sensitivity  calculation. 

element  mass  orig  =  element  mass;  %  Copy  properties  over  to  retain  them. 
element~EI_orig  =  elemental;  % 


%  Prompt  for  number  of  mode  frequencies  for  which  to  generate  sensitivities 

%  use  the  following  three  lines  for  user's  input  for  number  of  modal 
%  frequencies  for  which  to  generate  sensitivity  or  use  the  fourth  line 
%  which  uses  the  maxi  mum  number  of  modes  for  the  defined  system 

%d i  s p  ('max  n u mb e r  of  modes  '  ) 

%l  e  n  g  t  h  ( o  s  e  t ) 

%num  modesO  =  input)'  Enter  number  of  elastic  modes  for  sensitivity 
cal  cul  a t  T o n s  > >  '  ) ; 

num  modesO  =  ndof  ■  I  ength(BC(i  cnt  oset,:));  %ma  x  number  of  modes  in  system 
s  t  a  r  t  _  mo  d  e  =  num_rbm+  1;  %Skip"the  rigid  body  modes. 

%  Initializes  two  vectors  to  be  used  in  ploting  program. 
phiXPLOT  =  zeros(  ndof,  num  modesO);  %c  anti  I  ever  beam 
phi  AP LOT  =  z er  os  ( ndof ,  num" modes 0) ;  %c a n t  i  I  ever  beam 


MASS  SENSI T I  VI TY  CAL  CUL ATI  ON  LOOP 


sens  massO  =  0; 

if  massjbls  ~=  0;  %  from  Beam_XPrompt  as  user  inputs 

for  i  c  n  t  _  d  v  =  1:  num_el  ements;  %  loop  to  create  sensitivity  matrix 

%  Resetting  BeamX  properties  to  BeamA  properties, 
e  I  e  me  n  t  _  ma  s  s  ( :  ,  2 )  =  e  I  e  me  n  t  _  ma  s  s  (  :  ,  1 ) ; 

%  each  element,  one  at  a  time  will  have  a  change  in  mass 
element  ma  s  s  ( i  cnt  d  v ,  2 )  =  element  ma  s  s  ( i  cnt  d  v ,  2 )  *  ... 
(l+mass_change/100); 

A  s  s  e  mb  I  e  2  B  e  a  ms  _  c  r  s ;  %  Run  script  to  a  s  s  e  mb  I  e  b  e  a  ms . 

%  A  r  t i c i  a  I  Boundry  Conditions 
maO  base  =ma (BCoset,  BCoset); 

%  new  mass  matrix  of  the  system  defined  by  OSET 
mxO  =  mx(  BCoset ,  BCoset ) ; 

plotmx  =  mx_beam(  BCoset,  BCoset); 
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%  resulting  mass  matrix  with  ABC  used  in  plotting 
plotkx  =  kx  beam(BCoset,  BCoset); 

%  resul  ti  ng-sti  ffness  matrix  with  ABC  used  in  plotting 

k a 0_ b a s e  =  k a ( BCoset ,  BCoset ) ; 
kxO"=  k x ( BCoset ,  BCoset ) ; 

%  lam  (natural  freqA2,  r  a  d  A2/ s  ec  A2 ) ,  phi  (mode  shapes) 

[  I  amaOSET,  phi  aOSET]  =  f  Modes(  kaO  base,  maO  base) ; 

%  natural  freq  of  new  artifically  bounded  base  system 
[  I  amxOSET,  phi  xOSET]  =  f  Modes  ( k x 0,  mxO) ; 

%  natural  freq /  modes  of  new  artifically  bounded  system  with 
%  either  El  or  mass  changes  added  si  mi  liar  to  orginial  calculations 
[  I  amxpl  ot,  phi  xpl  ot]  =  fModes(  pi  otkx,  pi  otmx) ; 

%r  e  s  u  I  t  i  n  g  lam  &  phi  of  ABC  system  used  in  plotting 

%  Mode  shapes 

if  I  engt  h( pi  nned) ==1  &  pinned  ==  3  %  the  next  DOF  acts  a  little 

di  f  f  er  ent 

%  because  of  the  location  on  beam  thus  it  was  easier 
%  t  o  program  it  seperately 

phi  APLOT(  2:  ndof  -  2,  :)  =  phi aOSET( 1: ndof - 3,  :); 
phi XPLOT( 2:  ndof- 2,  :)  =  phi  xpl  ot ( 1:  ndof - 3,  :); 

elseif  I  engt  h(  pi  nned)  >1  &  pi  nned(  1,  1)  ==3 

phi  APLOT(  2:  pi  nned(  1,  2)  -  2,  :  )  =  phi  aOSET(  1:  pi  nned(  1,  2)-3,  :  ) 

phi  APLOT(  pi  nned(  1,  2)  -  1:  ndof  -  3,  :  )  =  phi  aOSET(pi  nned(l,  2)-2:  ndof-4,  :  ) 


elseif  length(pinned)<2 

phi APLOT( 1:  pi  nned- 3,  :)  =  phi  aOSET( 1:  pi  nned- 3,  :); 

phi  A P L OT ( pi  nned- 1:  ndof ■ 2, : )  =  phi  aOSETj pi  nned- 2 :  ndof - 3, : ) ; 

phi XPLOT( 1:  pi  nned- 3,  :)  =  phi  xpl  ot (  1:  pi  nned- 3,  :); 

phi  XPLOT(  pi  nned  -  1:  ndof  -  2,  :  )  =  phi  xpl  ot(pi  nned-2:  ndof-3,  :  ); 


else 

phi  APLOT( 1:  pi  nned(  1,  1)  -  3,  :)  =  phi  a  OS  ET(  1:  pi  nned(  1,  1)  -  3,  :); 

phi  APLOT(  pi  nned( 1, 1) -  1:  pi  nned(  1,  2)  -  3,  :  )  =  phi  a  OS  ET( pi  nned(  1,  1) - 

1:  pi  nned( 1,  2) -  3,  :  ) ; 

phi  APLOT(  pi  nned( 1,  2) -  1:  ndof  -  3,  :  )  =  phi  aOSET(  pi  nned( 1,  2) -  3:  ndof- 

4,:); 

phi  XPLOTf 1:  pi  nned(  1,  1)  ■  3,  :)  =  phi  xpl  ot ( 1:  pi  nned(  1,  1)  -  3,  :); 

phi  XPLOT(pi  nned(l,  1)-1:  pi  nned(l,  2)-3,  :  )  =  phi  xpl  ot(  pi  nned(  1,  1)- 

1:  pi  nned( 1,  2) - 3, :  )  ; 

phi  XPLOT( pi  nned( 1,  2) -  1:  ndof- 3, :  )  =  phi  xpl  ot ( pi  nned( 1,  2) - 3:  ndof - 

4,:); 

end  % 


num  rbmOSET  =  I  engt  h(  f  i  nd(  I  amaOSET  <  1)); 

%  f  T  n  d  number  of  rigid  bodies  in  new  ABC  system 
start_mode  =  num_ rbmOSET  +  1;  %  Skip  the  rigid  body  modes. 

faO  =  sqrt(  I  amaOSET)  /  ( 2*pi  ) ;  %n  a  t  u  r  a  I  freq  of  the  ABC  in  Hz 

%  F  o  r  m  ma  s  s  derivative  ma  t  r  i  c  e  s : 
m  deltaO  =  (mxO  -  maO  base)/(mass  change/100); 

%“ converts  percent  change  to  decimal  value 
%  NOTE:  mass  change  needs  to  be  the  same  change  used  in 
%  orginial  mass  sensitivity  matrix  calculation 

%  Mode  freq  sens  loop: 

end  mode  =  start  mode  +  (nummodesO  -  1); 

r  ow~num  =  0; 

for  icnt  modes  =  start  mode:end  mode; 
row  hum  =  row  num  +1; 

sens  massOfrow  num, icnt  dv)  =  phi  a  OS  E  T  (  :  ,  i  cnt  modes)1  * .  . 
[- 1  amaOSETCi  cnt  modes)  *  mdeltaO)  *... 
phi  a  OS  E  T ( : , i  c  nt "  modes ) ; 
end;  %  e  n  d  "for  i  c  n  t  _  mo  d  e  s 11  inner  loop 

end;  %  End  "for  i  c  n  t  _  d  v "  outer  loop  for  sensitivity  calculations 
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di  spl  acmentPI  ot_OSET  %  calls  a  program 

%  This  program  asse mb  les  the  beam  displace  me nt  vector  for 
%sens_beam(dispX)  and  base  beam(dispA)  under  ABC 

if  dispX  tot  ==  0;  %  when  the  vector  is  empty 

dispX"tot  =  displ;  %  displacement  vector  used  in  plotting 
di  spA"t  ot  =  di  spla; 

else 

dispX  tot  =  c  a  t  ( 1 ,  dispX  tot,  d  i  spl) ; 
dispA"tot  =  catjl,  dispA  tot,  di  spla) ; 
end  %  end"11  i  f  di  spX_t  ot  ==  0" 

end;  %  end  "if  ma  s  s  _  I  b  I  s  =  0 " 


0^  El  SENSI T I  VI TY  CALCULATI  ON  LOOP 

sens  El  0  =  0; 
i  f  Ef_ I  bl  s  ~=  0; 

for  i  c  n  t  _  d  v  =  1:  num_el  ements;  %  loop  to  create  sensitivity  matrix 

%  Resetting  BeamX  properties  to  BeamA  properties 
el  ement_EI  (  :  ,  2)  =  el  ement_EI  (:  , 

%  each  element,  one  at  a  time  will  have  a  change  in  El 
element  E I  ( i  c  n  t  d  v ,  2 )  =  element  E I  ( i  c  n  t  d  v ,  2 )  *  ... 

( 1  +  ( El  change/  100)  );' 

A  s  s  e  mb  I  e  2  B  e  a  ms  _  r  I  a ;  %  Run  script  to  a  s  s  e  mb  I  e  b  e  a  ms . 

%  A  r  t i f i  c  a  I  Boundary  Conditions 
maO  base  =  ma ( B C o s e t ,  BCoset); 
mxO"=  mx(  BCoset ,  BCoset ) ; 

plotmx  =  mx  beam(  BCoset,  BCoset); 
plotkx  =  kx“beam(  BCoset ,  BCoset); 

kaO  base  =  k a ( BCoset ,  BCoset ) ; 
kxO"=  k x ( BCoset ,  BCoset ) ; 

%  lam  (natural  freqA2,  r  a  d  A2/ s  ec  A2 ) ,  phi  (mode  shapes) 

[  I  amaOSET,  phi  aOSET]  =  f  Modes(  kaObase,  maO  base); 

[  I  amxOSET,  phi  xOSET]  =  f  Modes  ( k  x  07  mxO) ;  %sens  info 
[  I  amxpl  ot,  phi  xpl  ot]  =f  Mo  desjplotkx, plotmx);  %p  I  o  t  info 


if  I  engt  h(  pi  nned)  ==1  &  pinned  ==  3  %S  p  e  c  i  a  I  case  of  DOF  #3  pinned 
phi APLOT( 2:  ndof - 2,  :)  =  phi  a  OS  ET( 1:  ndof ■ 3,  :); 
phi XPLOT( 2:  ndof- 2,  :)  =  phi  xpl  ot ( 1:  ndof - 3,  :); 

elseif  I  engt  h(  pi  nned)  >1  &  pi  nned(  1,  1)  ==3 

phi  APLOT( 2:  pi  nned(  1,  2)  -  3,  :  )  =  phi  a  OS  ET(  2 :  pi  nned(  1,  2)  -  3,  :  ) ; 
phi  APLOT(  pi  nned( 1,  2) -  1:  ndof ■ 3,  :  )  =  phi  aOSET(  pi  nned( 1,  2) - 2:  ndof 


elseif  length(pinned)<2  %For  single  ABC 

phi  A  P  L  OT ( 1:  pi  nned- 3,  :)  =  phi  aOSET( 1:  pi  nned- 3,  :); 

phi APLOT( pi  nned- 1:  ndof- 2, : )  =  phi  aOSET( pi  nned- 2:  ndof- 3, : ) ; 

phi  XPL  OT  ( 1:  pi  nned-3,  :)  =  phi  xpl  ot  (  1:  pi  nned- 3,  :); 

phi XPLOT( pi  nned  - 1:  ndof- 2, : )  =  phi  xpl  ot ( pi  nned- 2 :  ndof - 3, : ) ; 

else  %F  o  r  t  wo  ABC's 

phi  APLOT( 1:  pi  nned(  1,  1)  -  3,  :)  =  phi  a  OS  ET( 1:  pi  nned(  1,  1)  -  3,  :); 

phi  APLOT(  pi  nned( 1, 1) -  1:  pi  nned(  1,  2)  -  3,  :  )  =  phi  a  OS  ET( pi  nned(  1,  1) 

1:  pi  nned( 1,  2) -  3,  :  ) ; 

phi  APLOT(  pi  nned( 1,  2) -  1:  ndof  -  3,  :  )  =  phi  aOSET(  pi  nned( 1,  2) -  2:  ndof 

4,:); 

phi  XPLOT( 1:  pi  nned(  1,  1)  -  3,  :)  =  phi  xpl  ot ( 1:  pi  nned(  1,  1)  -  3,  :); 

phi  XPLOT(pi  nned(l,  1)-1:  pi  nned(l,  2)-3,  :  )  =  phi  xpl  ot(  pi  nned(  1,  1) 

1:  pi  nned( 1,  2) - 3, :  )  ; 

phi  XPLOT( pi  nned( 1,  2) -  1:  ndof- 3, : )  =  phi  xpl  ot ( pi  nned( 1,  2) - 2:  ndof 

4,:); 

end 
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num  rbmOSET  =  I  engt  h(  f  i  nd(  I  amaOSET  <  1)); 

start_mode  =  num_ rbmOSET  +  1;  %  Skip  the  rigid  body  modes. 

faO  =  sqrt(  I  amaOSET)/(2*pi  ) ;  %n  a  t  u  r  a  I  freq  of  the  ABC,  Hz 
%  F  o  r  m  E I  d  e  r  i  v  a  t  i  v  e  ma  t  r  i  c  e  s : 

k_  d  e  I  t  a  0  =  (kxO  -  k  a  0_  b  a  s  e )  /  ( El  _  c  ha  nge/  100) ;  %  in  %/ 10  0 

%  Mode  freq  sens  loop: 

end  mode  =  start  mode  +  (nummodesO  -  1); 

r  o  w"  n  u  m  =  0; 

for  icnt  modes  =  start  mode:end  mode; 
row  iium  =  row  num  +1; 

sens  El  0(  row  iium, icnt  dv)  =  phi  aOSET( :  ,  i  cnt  modes)1  *... 
I  del  t  a  0” *  phi  a  OS  E  T ( : , i  cnt  modes); 
end;  %e  n  d " 11  f  o  r  icnt  modes" 


end;  %  End  "for  i  c  n  t  _  d  v 11  outer  loop  for  sensitivity  calculations 


di  spl  acmentPI  ot_OSET  %  cal  Is  a  program 

%  This  program  asse mb  les  the  beam  displace  me nt  vector  for 
%sens_beam(dispX)  and  base  beam(dispA)  under  ABC 

if  di  spX  t  ot  ==  0; 
di  s  pX"t  ot  =  di  s  pi; 
di  spA"t  ot  =  di  spla; 

else 

dispX  tot  =  cat(  1,  di  spX  tot,  d  i  spl) ; 
dispA"tot  =  c  a  t  ( 1 ,  dispA  tot,  di  spla); 
end  %  end"11  i  f  di  spX  tot  ==  []  " 

end;  %  e  n  d  "if  E I  _  I  b  I  s  =  [  ] " 


%  Copy  element  El  and  mass  properties  back  into  arrays: 
element  El  =  e  I  e  me  n  t  El  orig; 
el  ement“mass  =  element"mass_orig; 

clear  element  El  orig  element  mass  orig  end  mode  start  mode 
clear  row_num"icnt_dv  icnt_modes  Ibis  num_ET_dv  n  u  m_  ma  s  s  _  d  v 

%  asse  mb  I  e  of  total  sensitivity  matrix  for  ABC  System 
if  sens_massO  ==  0  &  s e n s _ E I  O  ~=0 ; 

T_  s  e  n  s  O  =  s  e  n  s  _  E I  O;  %  no  changes  in  mass 

elseif  sens  massO  ~=  0  &  sens  El  0  ==  0; 

T_sensO"=  sens_massO;  %  no  changes  in  El 


else 

T  sensO  =  c  a  t  ( 2 ,  sens  massO, sens  EIO); 

%" c h a n g e s  in  both  mass  and  El 

end  %end  "if  s e n s _ massO  ==  0  &  s e n s _ E I  0  ~=0 " 

%  Builds  complete  sens  matrix  of  all  ABC  systems 
if  T_sens_oset  ==  0;  %  when  mat  r  i  x  is  originally  empt  y 

T_sens_oset  =  T_  s  e  n  s  O; 

else  %  after  the  ma  t  r  i  x  has  s  o  me  values 

T_sens  oset  =  c a t ( 1 , T  sens  oset,  T  sensO); 
end  %end  "Tf  T_ s e n s _ os et " ==  [ I" 

I  amaOSET  =  I  a  ma  OS  E  T  ( f  i  nd(  I  amaOSET  >  1));  %s  kips  Rigid  body  modes 
v  e  c  t  _  I  a  mO  =  I  amaOSETf  1:  num_modesO) ;  %  nat  freq  of  base  beam  under  ABC 

%b u i  I  d s  a  vector  of  natural  freq  of  ABC  systems 
if  vect_lam_oset  =  =  0 ;  %  when  matrix  is  originally  empty 


vect  lam  oset  =  v e c t  lamO; 
e  I  s  e  %  after  the  ma  t  r  i  x  has  s  o  me  values 

vect  lam  oset  =  c a t ( 1 ,  vect  lam  oset,  vect  lamO); 
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end  %e  n  d  "if  v  e  c  t  lam  o  s  e  t  =  =  0 11 
dispC  ')  "  " 

oset  choice  =  input!1  Another  cycle  of  ASET  and  OSET  (  y  /  n ) ?  1  ,  '  s 1  )  ; 

%  loop  runs  until  11  n 11  is  inputed  creating  sens  matrix  &  lam  vector 
%  f  o  r  all  ABC  s  y  s  t  e  ms 
dispC  1  ) 

end;  %  end  "while  oset_choice  -='  n1  11 

%  ********************  BeamSensitivityOSET  rla  m  *********************** 
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BOUND ARYCONDITIONS  RLA 


BoundaryCondi  t i  ons_r I  a.  m 


%  Wr  i  1 1  en  by  Prof  Gor  di  s 

%  This  script  prompts  the  user  boundary  condition  information 
%  The  script  creates  a  vector  of  DOF  (with  respect  to  the  unrestrained 
%  structure)  and  then  extracts  the  rows  and  columns  of  the  complementary 
%  DOF. 

%  Script  defines  vector  "free  dof  set"  containing 
%  list  of  unrestrained  dof. 

%  The  boundary  conditions  are  applied  in  this  script. 

%  I  n  p  u  t  s  needed: 

% . 

%  ndof 

%  k a ,  ma,  kx;  mx 

%  Out  put  s : 

% . 

%  f  r  e  e  dof  set 
%  updated  Ka,  ma,  kx,  mx 
%  i  c  n  t  dof 
%  add  dof 
%  be  node 
%  b  c "  c  o  o  r  d 
%  be" DOF 
%  b  c  boolean 
%  alT  dof  s 
%  restrai  nt  swi  tch 

%  S  t  a  r  t  code: 


if  exi  st ( ' f r ee_dof  set 1 ) ==0;  %  Build  free_dof_set  vector 

d i  s p ( 1  Select  a  boundary  condition  set:1) 
disp('  (1)  Cl  amped-free1  ) 

dispj'  (2)  Cl  amped-CI  amped'  ) 

d  i  s  p  ( 1  (3)  Pinned-Pinned1) 

dispj1  (4)  User-  Def  i  n e d 1  ) 

dispj'  ( 5 )  Free-Free1) 

BCChoi  ce  =  i  n  p  u  t  ( 1  >>  Enter  choice:  '); 

if  BC  Choice  ==  1;  %  Cl  a  mped- f  r  ee 

free  dof  set  =  [3:ndof]; 
restrai  n t  _ s wi  tch  =  1  y1  ; 

elseif  BC_ Choice  ==  2;  %  Clamped-Clamped 

free  dof  set  =  [ 3 :  ndof -  2  ] ; 
restrai  nt_swi  tch  =  1  y 1  ; 

elseif  BC_Choi  ce  ==  3;  %  Pinned-Pinned 

free  dof  set  =  [ 2 :  ndof -  2  ndof ] ; 
restrai  nt_swi  tch  =  1  y1  ; 

elseif  BC_Choi  ce  ==  4;  %  User-Defined 

i  c nt  dof  =0; 

add  dof  =  1  y 1 ; 

whi  Te  add  dof  ==  1  y 1 ; 


b c _ n o d e  =  input!1  Node  number  for  restraint  ?  "0"  to  end:  '); 


if  be  node  ==  0; 
break 

end; 

bc_coord  =  input!1  Translation  or  Rotation  ?  (t/r)  ','s'); 

i  c nt  dof  =  i  c nt  dof  +1; 
if  be  coord  ==  rt 1  ; 

be  DOF ( i  c nt  dof )  =2  *  be  node  -  1; 
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el  sei  f  be  coord  ==  1  r 1 ; 

be  DOF ( i  c nt  dof )  =  2  *  be  node; 
end;  ”  %  End  i  f :  e  I  s  e  block 

end;  %  End  while  a d d _ d o f 

be  boolean  =  ones(ndof,l);  %  [  1  1  1  .  .  .  ient  dof] 

bc"bool  e  a  n ( be  DOF)  =  z er os ( I  engt h(  be  DOF ) ,  1 ) ;  %  Put  zeros  in  restrained  dof 
alf  dofs  =  [l“ndof];  "  %List  of  all  dof 

free  dof  set  =  all  dofs(  I  ogi  cal  ( be  boolean));  %  Extract  free  dof 
r  est  r  ai  n t  _ s wi  t  c  h  =- 1  y 1  ; 

el  seif  BC_  Choice  ==  5;  %  Free  - free  beam 

free  dof  set  =  [  1 :  ndof  ] ; 
restrai  nt_swi  tch  =  1  n1  ; 

end;  %  End  if-elseif  choice  block 

end;  %  End  exist  block 

ka  =  kaffree  dof  set, free  dof  set); 
ma  =  ma(free"dof"set,free"dof"set); 
kx  =  k x ( f r ee"dof "s et , f r ee"dof "s et ) ; 
mx  =  mx(free~dof~set,free“dof~set); 

0/0*******>m*******end  BoundaryCondi  ti  ons  rla.m  *********************** 
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H.  BUILD2BEAMS  RLA 


B  u  i  I  d  2  B  e  a  ms  rla.m 


clear 
c  I  c 

%  R  e  v  i  s  i  o  n  history: 

% . 


%  Ver. 

1.  0: 

9/22/  94 

Basic 

t  wo  beam  ass embl  y 

% 

2.  0: 

Added 

mu  1  t  i  -  e  1  e  me  n  t  changes 

% 

2.  1 

3/28/  95 

Added 

read/write  to  file,  rebuild 

capability 

% 

2.  2 

3/29/  95 

Added 

1  u mp e d  mass  additions 

% 

3/ 10/ 04 

Added 

Sensitivity  ma  trices,  error 

pr  edi  ct i  on 

plots 


%  Program  Des  c  r  i  pt  i  o  n : 


%  This  program  asse mb les  the  mass  and  stiffness  ma trices  for  2  free-free 
%  b  e  a  ms ,  referred  to  as  11 B  e  a  mA 11  (analysis)  and  11 B  e  a  mX 11  ( e  x  p  e  r  i  me  n  t  a  I  ) .  The 
%  program  can  be  run  in  several  modes: 

% 

%  "Build"  mo  d  e : 


%  The  user  provides  baseline  data  for  BeamA,  assumed  to  be  a 
%  homogeneous,  uniform  beam.  Data  provided: 

% 

%  (1)  Beamlength 

%  ( 2 )  N u mb e r  of  e I  e me n t  s 

%  (3)  Nominal  El 

%  (  4 )  Nominal  cross-sectional  area 

%  (5)  No  mi  nal  wei  ght  density 

% 

%  The  program  then  prompts  the  user  for  instructions  on  howto  modify 
%  "  B  e  a  mA "  data  to  arrive  at  "  B  e  a  mX 11  data.  The  user  can  mo  d  i  f  y  e  I  e  me  n  t 
%  masses,  and/or  element  El  values.  The  modification  can  be  applied  to 
%  either  a  single  element,  or  range  of  elements,  e.g. 

% 

%  Modify  single/range  element  mass  values  (  y  /  n )  ?  y 

% 

%  If  11  y "  is  entered,  the  user  enters  the  n u mb e r  of  the  e I  e me n t  for  ma s s 
%  a  d  j  u  s  t  me  n  t : 

% 

%  Enter  e I  e me n t  I  a b e I  ( s )  for  mass  mo d i  f  i  c a t  i  o n :  1 

%  Use  MAT  LAB  vector  f  o  r  ma  t  >  1  3  5:7  9 

% 

%  Enter  percentage  mass  change  ( +/ -  %) 

% 

%  The  user  is  pro mp ted  to  modify  another  element  or  range  of  elements: 

% 

%  Modify  another  element  mass  value  (  y  /  n )  ?  y 

% 

%  This  process  continues  until  the  user  enters  an  11  n 11  for  no  change. 

%  This  entire  process  can  then  be  repeated  for  El  adjust  me  nt. 

% 

%  The  program  saves  the  beam  definition  data  in  a  binary  (.mat)  file 
%  "beamdata"  at  the  end  of  execution. 

% 

%  The  program  can  also  be  run  in  "Read"  mode  by  entering  an  "r"  at 
%  the  initial  p  r  o  mp  t . 

% 

% 

%  Script  Execution  Path: 

% . 


%  B  u  i  I  d  2  B  e  a  ms  c  r  s .  m 

%  B  e  a  mA  P  r  o  mp  t  c  r  s .  m 

%  BeamX- Pr  ompt"cr  s.  m 

%  As  s  e  mb  I  e2  Be  a  ms  c  r  s .  m 

%  "  [  mx ] 

%  AddLumpmass  crs.  m 

%  BoundaryConditions  crs.m  -- 

%  Plot  Be  a  mModes  crs.m 

% 


User  executes  this  program. 

Prompts  User  for  BeamA  nominal  beamdata 
Prompts  User  for  BeamX  modification  beamdata 
Called  by  Bui  I  d2Beams,  builds  [ka]  [ma]  [kx] 
plots  f  r  eqs . 

--  Prompts  User  for  BeamX  lumped  mass  addition 
Prompts  user  for  B.C.'s  and  applies  them. 
Calculate  beam  modes  and  plot  frequencies 
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%  BeamSensitivity  crs.m  --  Calculate  sensitivity  matrix  T-sens 

%  Beams ensitivityOSET  crs.m--  Calculate  sensitivity  matrix  using  ABC 

%  recorded  H  crs.m  "  --  Calculates  the  nat.  freq  of  BeamX  with  ABC 

applied 

%  Assemble  Sens  crs.m  --  Assembles  the  sens  matrices  and  calculates 

errors. 

%  ABCrunTHRU.m  --  Calculates  the  DV  and  cond  number  of  matrix  used 

%  Saves  data  to  11  b  e  a  md  a  t  a .  ma  t 11 


% 

%  Start  code: 


d i s  p ( '  Building  2  beams  from  scratch...1) 

BeamA  Prompt  rla;  %  Prompt  for  BeamA  Data:  run  prompt  script 

Be  a  mX"  Pr  ompt  "  r  I  a ;  %  Prompt  for  BeamX  Modification  Data: 

A  s  s  e  mb  I  e  2  B  e  a  ms  rla;  %  Run  script  to  assemble  mass  and  stiffness  ma  t  r  i  c  e  s 

A  d  d  L  u  mp  ma  s  s  rla;  %  B  e  a  mX  I  u  mp  e  d  mass  vector  construction  and 

%  "  application 

kx  beam  =  kx;  %  saves  the  BeamX  matrices  without  BC  to  be  used  later 
mx“beam  =  mx; 

BoundaryCondi ti  ons_rl  a;  %  Prompt  for,  and  apply  boundary  conditions 

kx  beamBC  =  kx;  %  saves  the  BeamX  matrices  with  BC  to  be  used  later 
mx“beamBC  =  mx; 

ka“beamBC  =  ka;  %  saves  the  Beam  A  matrices  with  BC  to  be  used  later 
ma“  beamBC  =  ma; 

PI  o  t  B  e  a  mMo  d  e  s  _  r  I  a  %  Calculate  beam  modes  and  plot  frequencies 

BeamSensitivity  rla;  %  Calculate  sensitivity  matrix  T-sens 
BeamSensitivity  OS  ET_  rla;  %  Calculate  sensitivity  matrix  using  ABC 


recorded  H  rla;  %  C  a  I  u  I  a  t  e  s  the  nat.  freq  of  BeamX  with  ABC  applied 

Assembl  eSens  rla;  %  Assembles  the  sens  matrices  and  calculates  errors. 
ABCr  unTHRU  r  F  a ;  %  Calculates  the  DV  and  cond  number  of  matrix  used 

FOM  rla;  "  %C  a  I  c  u  I  a  t  e  s  the  Figure  of  Merit  for  each  prediction 


%  Save  Defining  Parameters  for  Beams  and  plots 

d i s p ( 1  ...saving  beam  data  to  file1) 
save  b  e  a  md  a  t  a .  ma  t 


d i  s  p ( 1  Bui  I  d2Beams  end.1) 

% 


END  Bui  I  d2Beams  rla.m 
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DAMAGE  PREDICTOR 


%  %  Mode  error  predictor 
% 


f  o  r  ma  t  long 

%///////////  ORGI  NAL  METHOD/ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
%c  a  I  c  u  I  a  t  e  s  d  v  for  base  system  s  u  mme  d 
for  q  =1 : 2  0 

dv  calculated  b  a  s  e  (  :  ,  q )  =T  sens  tot(  1:  q,  :  )\vect  lam  tot(l:q); 
end 

%c  a  I  c  u  I  a  t  e  s  dv  for  ABC1  system  su  mme  d 
for  t  =40: 59 

dv  calculated  basel(:,(t-39))  =  T_sens  tot(40:t,:)\vect  lam  t  ot  (  40:  t ) ; 
end 

%c  a  I  c  u  I  a  t  e  s  dv  for  ABC2  system  su  mme  d 
for  t  =60: 79 

dv  calculated  base2(:,(t-59))  =  T  sens  tot(60:t,:)\vect  lam  tot(60:t); 
end 

%c  a  I  c  u  I  a  t  e  s  dv  for  ABC3  system  su  mme  d 
for  t  =80: 99 

dv  calculated  bas  e  3  ( :  ,  ( t  -  79) )  =T  sens  t  ot  ( 80:  t ,  :  )  \  vect  I  a  m  t  ot  ( 80 :  t ) ; 
end 

for  t  =100: 119 

dv  calculated  bas e  3 ( : , ( t - 99) )  =T  sens  t ot (  100: t ,:) \ vect  I  a m  t ot (  100: t ) ; 
end 

%F  o  r  ABC  System  1  and  2  %%%%%%%%%%%%%%%%%%%%%%0/o%%%y  o  lo  c  omen  tee 


%  %%%%%%%%%%%%%  2  ABC  system  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%% 

%  T  c  o  mb  o  1  =  z  e  r  o  s  ( 9 ,  2  0 ) ; 

%  T~combol(l:3,:)  =  T  sens  t  o  t  (  1 :  3 ,  :  ) ; 

%  T“ c o mb o 1  (  4:  6,  : ) =  T“s ens't  ot ( 40: 42, : ) ; 

%  T” c o mbo 1 ( 7:  9,  : ) =  T"sens"t ot (  7  8: 8  0, : ) ; 

%  "  " 

%  vect  c  o  mb  o 1  =  z  e  r  o  s ( 9 ,  1 ) ; 

%vect"co mb 01(1:3,1)  =  vect  lam  tot(l:3,l); 

%  vect~combol(4:6,l)  =  vect~lam'tot(40:42,l); 

%  v  e  c  t  ~  c  o  mb  o  1  ( 7 :  9 ,  1 )  =  vect~lam'tot(78:80,l); 

% 

%  f  o  r  e  =  1 :  9 

%  p  c  o  mb  o  1  ( e ,  :  )  =  T  combol(l:e,:)\vect  combol(l:e,l) 

%  end 

%  b  b  b  b  =  ( p  c  o  mb  o  1 )  1 

%  %%%%%%%%%%%%  3  ABC  system  %%%%%%%%%%%%%%%%%%0/o0/o0/o%%% 

%  T  c o mb o 1  =  zeros(  12, 20) ; 

%  T” c o mb o 1 (  1:  3,  :  )  =  T  sens  t  ot ( 1:  3, :  ) ; 

%  T~ c o mbo 1 ( 4 :  6, :  )  =  T"s ens't  ot ( 40: 42, : ) ; 

%  T"combol(  7:  9,  :  )  =  T's ens't  ot  (  7  8:  8  0,  :  )  ; 

%  T'combolt  10:  12,  :  )  ="T  sens  t  ot  (  117:  119,  :  )  ; 

%  "  "  " 

%  vect  c  o  mb  o 1  =  z  e  r  o  s (  1 2 ,  1 ) ; 

%vect"co mb 01(1:3,1)  =  vect  lam  tot(l:3,l); 

%  vect"combol(4:6,l)  =  vect'l  am'  tot(40:  42,  1) ; 

%  vect“combol(7:  9,  1)  =vect"lam'tot(78:80,l); 

%  v  e  c  t "  c  o  mb  o  1  ( 10:  12,  1 )  =  vect  lam  tot(117:119,l); 

% 

%  f  o  r  e  =  1:12 

%  p  c  o  mb  o  1  ( e ,  :  )  =  T_combol(l:e,:)\vect  combol(  1:  e,  1) 

%  end 
% 

%  b  b  b  b  =  ( p  c  o  mb  o  1 ) 1 

%%%%%%%%%%%  5  ABC  system  %%%%%%%%%%%%%%%%%%%%%%0/o0/o%%% 

T  combol  =  zeros(  18,  20) ; 

T “ c o mb o  1  (  1:3,:)=  T  sens  t  o t  (  1 :  3 ,  :  ) ; 

T'combol(4:6,:)  =  T"sens~tot(40:42,:); 

T"combol(7:9,:)  =  T~sens~tot(78:80,:j; 

T'comboll  10:  12,  :  )  ="T  sens  t  ot  (  117:  119,  :  ) ; 

T'comboK  13:  15,  :  )  =  T's  ens't  ot  (  1  5  6:  1  5  8,  :) ; 
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T_combol(  16: 18, : ) =  T_sens_t  ot (  1  9  5:  1  9  7, :  ) ; 

vect  combol  =  zeros(  18,  1) ; 
vect"co mb 01(1:3,1)  =  vect  lam  tot(l:3,l); 
vect~combol(4:6,l)  =  vect~lam~tot(40:42,l); 
vect“co  mb  01(7:9,1)  =  vect~lam~tot(78:80,l); 
vect“combol(  10:  12,  1)  =  vect  lam  tot(117:119,l); 
vect"co mb o  1  ( 13:  15,  1 )  =  vect-l  am_tot(  156:  158,  1) ; 
vect^co mb ol(16:18,l)  =  vect“l  am“tot(  195:  197,  1) ; 

for  e  =1:18 

p  c  o  mb  o  1  (  e ,  :  )  =T  combol(l:e,:)\vect  combol(  1:  e,  1) 
end 

bbbb  =  ( p c o mb o 1 ) 1 


%  %For  base  and  ABC  System  %%%%%%%%%%%%%%%%%%%%%0/o0/o%%%y o  lo 
% 


%  %A  B  C  1  system 
%  T  c  o  mb  o  1  =  z  e  r  o  s  ( 6 ,  2  0 ) ; 

%  T~combol(l:3,:)  =  T  sens  t  ot  (  1:  3,  :  ) ; 

%  T“ c o mb o 1  (  4:  6,  : ) =  T“s ens“t  ot ( 40: 42, : ) ; 

%  "  " 

%  vect  c  o  mb  o 1  =  z  e  r  o  s ( 6 , 1 ) ; 

%vect"co mb ol(l:3,l)  =  vect  lam  tot(l:3,l); 

%  vect"combol(4:6,l)  =  vect"lam"tot(40:42,l) 

% 

%  f  o  r  e  =  1:6 

%  p  c  o  mb  o  1  ( e ,  :  )  =  T  combol(l:e,:)\vect  combol(  1:  e,  1) 
%  end 
% 

%  bbbb  =  ( p  c  o  mb  o  1 ) 1 
% 

%  %A  B  C  2  system 
%  T  combo2  =zeros(6,20); 

%  T“ c o mb o 2 (  1:  3,  :  )  =  T  sens  t  ot ( 1:  3, :  ) ; 

%  T “ c o mb o 2 (  4 :  6 ,  :  )  =  T“sens“t  ot (  7  8:  8  0,  :  )  ; 

%  "  " 

%  vect  c  o  mb  o  2  =  z  e  r  o  s ( 6 , 1 ) ; 

%vect"co mb 02(1:3,1)  =  vect  lam  tot(l:3,l); 

%  vect"combo2(4:6,l)  =  vect"lam"tot(78:80,l) 

% 

%  f  o  r  e  =  1:6 

%  p  c  o  mb  o  2  ( e ,  :  )  =T  combo2(l:e,:)\vect  combo2(  1:  e,  1) 
%  end 

%  c  c  c  c  c  =  c  o  n  d  (  T  c  o  mb  o ) 

%  c  c  c  c  =  ( p  c  o  mb  o  2 ) 1 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%y o  I  o  c  o  me  n  t  e  e 


c  o  me  n  t  e  e 
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J. 


DISPLACEMENTPLOT  OSET  RLA 


%  ********************  displace  me ntPlot  OSET  rla.m  *********************** 

%  Wr  i  1 1  e  n  by  Constance  Fernandez  Spring  2004 
%  Updated  by  John  Mentzer,  Spring  2007 

%  This  program  plots  the  mode  shapes  (phi,  lam)  phi  vs  nodal  position 
%  of  beam  when  ABC  are  applied. 

%  Inputs 

% . 

%  p  I  o  t  k  x 
%  kaO  base 
%  n um"  modes  0 
%  phi  XPLOT 
%  phi APLOT 
%  pinned 
%  num_el  ements 

%  Out  put  s 

% . 

%d  i  s  p  1 ,  d  i  s  p  I  a 

%j  j .  g 

%  ypos 
% . 

%  displ  =  zerosfcei  I  (.  5*si  ze(pl  otkx,  1)),  num  modesO); 

%i  n  i  t  i  I  i  z  e  disp  vector  and  provides  the  first  zero  of  the  vector, 
displ  =  zeros(cei  I  ((20/38)*si  ze(pl  otkx,  1)),  num  modesO); 

%i  n  i  t  i  I  i  z  e  disp  vector  and  provides  the  first  zero  of  the  vector. 

%  displa  =  zeros(ceil(.5*size(kaO  base, 1)), num  modesO); 

%i  n  i  t  i  I  i  z  e  disp  vector  and  provides  the  first  zero  of  the  vector, 
displa  =  zer  os  ( cei  I  ( (  2  0/  3  8  )  *s  i  ze(  kaO  base, 1)), num  modesO); 

%i  n  i  t  i  I  i  z  e  disp  vector  and  provides  the  first  zero  of  the  vector, 
for  jj  =  1:  cei  I  ((20/38)*si  ze(pl  otkx,  1))  ; 

d  i  s  p  1  ( j  j  +1 ,  :  )  =  phi  XPLOT(2*jj  -1,  1:  num  modesO); 

%  every  other  phi  to  give  displace  me  n  t  at  sequential  nodes 
di  s  p  1  a  ( j  j  +1,  :  )  =  phi  APLOT(2*jj  -1,  1:  num  modesO); 
end 

%  This  loop  normalizes  the  modes  shapes  to  the  tip  modal  displace  me  nt. 
if  pinned  ==41  %  t  i  p  pinned  is  a  special  case,  no  new  calculations  are  needed 
di  s  p  1  ( :  ,  : )  =  di  s  p 1 ( :  ,  :  ) ; 
di  s  p 1  a ( :  ,  : )  =  di  s  p 1  a ( :  ,  : ) ; 

elseif  I  ength(  pi  nned)  >1  &  pi  nned(  1,  2)  ==41  %t  i  p  pinned  with  2  ABC's 
di  s  p  1  ( :  ,  :  )  =  di  s  p  1  ( :  ,  :  ) ; 
di  s  p 1  a ( :  ,  :  )  =  di  s  p 1  a ( :  ,  :  ) ; 
else 

for  g  =  1 :  n u m  modesO 

displ(:,g)~=  displ(:,g)/displ(num  el  ements+1,  g) ; 
d  i  s  p  1  a  ( :  ,  g )  =  di  spla(  :  ,  g)/di  spla(num  el  ements+1,  g) ; 
end 
end 
%  end 

%y  p  o  s  =  [  1 :  1 :  n  u  m  el  ements+1] ;  %  Location  of  nodes  used  in  plotting 
ypos  =  [  1:  1:  num  el  ements+1];  %  Location  of  nodes  used  in  plotting 
%  i  f  mass  Ibis  ~=  [  ] ; 

% 

%  for  k  k  =  1 :  s  i  z  e  (  ma  s  s  I  b  I  s ,  1 ) ; 

%  ff  =0; 

%  for  JJ  =  1: I  engt h( f i  nd(  mass  I  bl  s  ( kk,  :  )  >0) ) ; 

%  ff  =  ff+1; 

%  p o s  m(  k k ,  2* J  J  - 1 )  =  ma  s  s  I  b I  s  ( k  k ,  f  f ) ; 

%  posm(kk,  2  *  j  j  )  =  mass  I  bl  s  ( kk,  f  f ) +1; 

%  end 

%  end 

% 

%  i f  kk  ==  1 
%  p  o  s  M  =  p  o  s  m; 

% 

%  else 

% 

%  f  o  r  u  u  =  1 :  k  k  -  1 ; 

%  posM  =  cat(2,  p  o  s  m(  u  u ,  :  ) ,  posm(  uu+1,  :  ) ) ; 

%  end 

% 

%  end 

%  posM  =  sort  (  p  o  s  M(  f  i  nd(  posM>0)  ) )  ; 
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%  m  =  .  5*ones(  si  ze(  posM) )  +i  cnt  oset; 

%  e  n  d 
% 

%  i  f  El  Ibis  ~=  [  ] ; 

%  for  kk  =  1 :  si  z  e  ( El  lbls,l); 

%  ff  =0; 

%  for  J  J  =l:length(find(EI  I  bl  s  ( kk,  :  )  >0) ) ; 

%  ff  =  ff+1; 

%  p o s e ( kk,  2*j  j  - 1)  =  El  I  bl  s( kk,  f f ) ; 

%  pos e(  kk,  2*j  j  )  =  El  I  El  s( kk,  f f )  +1; 

%  end 

%  end 

% 

%  | f  kk  ==  1 

%  p  o  s  E  =  p  o  s  e ; 

% 

%  else 

% 

%  f  o  r  u  u  =  1 :  k  k  -  1 ; 

%  posE  =  cat(2,  p  o  s  e  ( u  u ,  :  ) ,  pose( uu+1, : ) ) ; 

%  end 

%  end 

%  posE  =  sort(  posEffi  nd( posE>0) ) ) ; 

%  e  =  - .  5*ones(  si  ze( posE) ) +i  cnt  oset; 

%  end 


o/0  ********************  end  di  spl  acementPI  ot  OSET  rla.m  *********************** 
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K 


FBEAMKM  RLA 


%  %  ********************  fbeamkm  rla.m  *********************** 

% 

%  %  function  [  k  b  e  a  m,  mb  e  a  m]  =f  b  e  a  mk  m(  I  ,  e  i  ,  m) 

%  %  Pr  ovi  ded  by  Prof  Gor  di  s 
% 

%  function  [  k  b  e  a  m,  mb  e  a  m]  =f  b  e  a  mk  m(  I  ,  e  i  ,  m) 

%  % 

%  % 

%  %  This  function  returns  the  stiffness  and  mass  ma trices  for 
%  %  a  s  i  mp  I  e  2  -  n  o  d  e  b  e  a  m  e  I  e  me  n  t . 

%  % 

%  %  N  o  t  e :  m  =  r  h  o  *  area  *  length  =  t  o  t  a  I  e  I  e  me  n  t  mass 

%  % 

%  %  Reference:  R.D.  Cook,  Concepts  and  Applications  of  F.E.  Analysis 
% 

%  %  Out  put  s 
%  % . 

%  %  k  b  e  a  m,  mb  e  a  m,  i  ,  j 
% 

% 

%  kbea  m=z er  os ( 4,  4) ; 

%  mbea  m=z  er  os  (  4,  4) ; 

%  % 

%  kbea m(  1,  1)  =12.  0; 

%  k b e a m(  1,  2)  =6.  0*1  ; 

%  k b e a m(  1,  3)  =- 12.  0; 

%  kbea m(  1,  4)  =6.  0*1  ; 

%  k b e a m(  2,  2)  =4.  0*1  A2; 

%  k b e a m(  2,  3)  =-  6.  0*1  ; 

%  kbea m(  2,  4)  =2.  0*1  A2; 

%  k b e a m(  3,  3)  =12.  0; 

%  k b e a m(  3,  4)  =-  6.  0*1  ; 

%  kbea m(  4,  4)  =4.  0*1  A2; 

%  % 

%  mbea m(  1,  1)  =156.  0; 

%  mbea m(  1,  2)  =22.  0*1  ; 

%  mbea m(  1,  3)  =5  4.  0; 

%  mbea m(  1,  4)  =-13.  0*1  ; 

%  mbea m(  2 ,  2)  =4.  0*1  A2; 

%  mb e a m(  2 ,  3)  =13.  0*1  ; 

%  mbea m(  2 ,  4)  =-  3.  0*1  A2; 

%  mbea m(  3,  3)  =156.  0; 

%  mb e a m(  3,  4)  =-22.  0*1  ; 

%  mbea m(  4,  4)  =4.  0*1  A2; 

%  % 

%  for  i  =1:  4; 

%  for  j  =i  :  4; 

%  kbea m(  j  ,  i  )  =k b e a  m(  i  ,  j  ) ; 

%  mbea m(  j  ,  i  )  =mbea m(  i  ,  j  ) ; 

%  end 
%  end 
%  % 

%  k  b  e  a  m=(  ei  / 1  A3)  *  k  b e  a  m; 

%  mb  e  a  m=(  m/  4  2  0 .  0  )  *  mb  e  a  m; 

%  % 

%  %  end  function  b  e  a  mk  m 
% 

%  %  ********************  END  fbeamkm  rla.m  *********************** 
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L. 


FMODES 


%  /o  f  M°deS' m 
%  function  [  I  a m,  ph i  ]  =f  modes  ( k,  m,  numto  print); 

%  %  Pr  ovi  ded  by  Prof  Gor  di  s 

%  %  This  program  prints  to  the  screen  natural  modes  of  system  (phi). 

%  % 

%  %  Usage:  [I  am,  phi  ]=fmodes(k,  m,  num  to  print) 

%  %  "  " 

%%  This  function  can  be  used  with  1  to  3  arguments,  as  follows: 

%  % 

%%  [  I  am,  phi  ]  =f  modes  ( a)  :  Produces  modes  of  [a]  with  no  print  of  freqs  in 

Hz. 

%%  [  I  am,  phi  ]  =f  modes  ( a,  i  )  :  Produces  modes  of  [a]  with  print  of  11  i  11  freqs  in 

Hz . 

%%  [  I  am,  phi  ]  =f  modes  ( k,  m)  :  Produces  modes  of  [m\k]  with  no  print  of  freqs  in 

Hz. 

%  % 

%  % 

%  % 

%%  This  function  returns  a  vector  containing  eigenvalues  (rad/ sec)  A2 

%  %  and  a  ma t  r  i  x  containing  the  mass  n o r  ma I  i  z e d  mode  shapes. 

%  %  The  mode  information  is  sorted  by  frequency  in  ascending  order. 

%  %  If  numto  print  >0;  tabular  listing  of  numto  print  freqs  in  Hz  is 

printed. 

%  %  If  numto  print  <=0,  no  print. 

%  "  " 

%  %  Inputs 

%  % . 

%  %  v,  index,  m,  k 
% 

%  %  P  r  o  g  r  a  ms 

%  % . 

%  %  f  N  o  r  ma  I  i  z  e 
% 

%  %  Out  put  s 

%  % . 

%  %  phi 

%  %  num  t  o  _  p  r  i  n  t 
%  %  e  r  r  o  r 
%  %  e 
% 

%  % 


% 

%  i  f  n  a  r  g i  n  ==  1 ; 

%  [A]  w/  no  print  request  for  freqs  in  Hz. 

% 

%  v  ( 1 ,  :  )  =  1  n  o  r  ma  I  i  z  a  t  i  o  n 
%  [  v,  d]  =ei  g(  k) ; 

%  [temp,  i  ndi  ces]  =  s  o  r  t  ( a  b  s  ( d  i  a  g  ( d ) ) ) ; 

%  lam  =  diag(d); 

%  lam  =  lam(indices); 

%  [phi]=fNormalize(v(:, indices),  'one'); 

%  num  to  print  =  0; 

%  "  " 

%  e  I  s  e  i  f  n  a  r  g  i  n  ==  2  &  s  i  z  e  (  m,  1 )  ==  1 ;  %  [A]  w/  print 

request  for  freqs  in  Hz. 

% 

%  v  ( 1 ,  :  )  =  1  n  o  r  ma  I  i  z  a  t  i  o  n 
%  [  v,  d]  =ei  g(  k) ; 

%  [temp,  i  ndi  ces]  =  s  o  r  t  ( a  b  s  ( d  i  a  g  ( d ) ) ) ; 

%  I  a  m  =  d  i  a  g  ( d ) ; 

%  lam  =  lam(indices); 

%  [phi]=fNormalize(v(:, indices),  'one'); 

%  n  u  m  t  o  p  r  i  n  t  =  m; 

%  "  " 

%  e  I  s  e  i  f  n  a  r  g  i  n  ==  2  &  s  i  z  e  (  m,  1 )  >  1 ;  %  [  k  ] ,  [  m] 

w/  no  print  request  for  freqs  in  Hz. 

% 

%  mass  n o r  ma  I  i  z a t  i  o n 
%  [  v,  d]  =ei  g(  m\  k) ; 

%  [lam,  index]  =sort(abs(diag(d))); 

%  [phi  ]=fNormal  i  ze(v(:  ,  i  ndex),  '  mass'  ,  m); 

%  n  u  m  t  o  p  r  i  n  t  =  0 ; 

%  ~  " 
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[  k  ] ,  [  m]  w/  print 


%  e  I  s  e  i  f  n  a  r  g  i  n  ==  3  &  s  i  z  e  ( k ,  1 )  >  1  &  s  i  z  e  (  m,  1 )  >  1 ;  % 
request  for  freqs  in  Hz. 

% 

%  mass  n o r  ma I  i  z a t  i  o n 
%  [  v,  d]  =ei  g(  m\  k) ; 

%  [lam,  index]  =sort(abs(diag(d))); 

%  [phi  ]=fNormal  i  ze(v(:  f  i  ndex),  1  mass'  ,  m); 

% 

%  else 
% 

%  n  u  m  t  o  p  r  i  n  t  =  - 1 ; 

%  error(rError  from  fModes.m:  Check  input  arguments.1) 

% 

%  end 
% 

%if  numto  print  >  I  e  n  g  t  h  ( k ) ; 

%  "numto  print  =  I  engt  h(  k) ; 

%  e  n  d 
% 

%  if  nargin  <  3  &  rem(l  ength(k),  2)  ==0  &  k(l:  I  ength(k)/2,  1:  I  ength(k)/2)  == 
zeros(length(k)/2,length(k)/2);  %Have  [A]  matrix 
%  e  =  1;  %  Ei  genval  ues  ar  e  wn 

%  else 

%  e  =  0.5;  %  Eigenvalues  are  wn^2 

%  end 
% 

% 

% 

%if  num  to  pri  nt  >0; 

% 

%  d  i  s  p  ( 1  1  ) ,  d  i  s  p  ( 1  1  ) 

%  dispC . ') 

%  disp('Freqs  in  Hz.:1) 

%  disp((lam(l:num  to  print).  Ae )  /  2  /  p  i  ) 

%  dispC . . ) 

% 

%  end 
% 

%  ********************  END  fModes.m  *********************** 
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M.  FNORMALIZE 


f  N  o  r  ma I  i  z  e .  m 


%  function  [phi]  =  fNormal  i  ze(phi  f  method,  m); 
%  % 


%  %  Usage:  [phi]  =  fNormal  i  ze(  phi  ,  method,  m) ; 

%  % 

%%  phi:  matrix  whose  columns  are  to  be  (independently)  normalized 
%%  method:  String  variable.  The  following  choices  are  available: 

%  % 


%  % 


%  % 
%  % 


mass' 
inf1 
one1 
I  e  n  g  t  h 1 


Mass  nor  mal  i  zat i  on 
Infinity  n  o  r  ma  I  i  z  a  t  i  o  n 
First  element  =1 
Length  =  1 


%  %  m:  ma  t  r  i  x  used  for  normalization,  i  .  e .  ,  phi  1  *  m*  p  h  i  =  eye 
%  % 


s  wi  t  c h  met  hod 

case  'mass' 


Mass  nor  mal  i  zat  i  on 


%  %  d  i  s  p  ( '  ma  s  s  n  o  r  ma  I  i  z  a  t  i  o  n 1  ) 

%  phi  =  phi  *diag(sqrt(diag((phi'  *  m  *  p  h  i  ) .  A(  - 

% 

%  c  a  s  e  '  i  n  f '  %  I  n  f  i  n  i  t  y  n  o  r  ma  I  i  z  a  t  i  o  n 

%  %  disp('  i  nf  normal  ization1  ) 

%  f  o  r  i  c  n  t  c  o  I  s  =  1 :  s  i  z  e  ( p  h  i  ,  2 ) ; 

%  pi  i  ( :  ,  i  c  n  t  c  o  I  s )  = 

phi(:,icnt  cols)/norm(phi(:,icnt  cols),  inf); 

%  "  end 

% 

%  c  a  s  e  '  o  n  e '  %  F  i  r  s  t  e  I  e  me  n  t  =  1 

%  %  d  i  s  p  ( '  o  n  e  n  o  r  ma  I  i  z  a  t  i  o  n 1  ) 

%  phi  =  phi  *  di  a g ( ( phi  ( 1,  :  ) .  A( - 1) )  '  ) ; 

% 

%  c  a  s  e  '  I  e  n  g  t  h '  %  L  e  n  g  t  h  =  1 

%  %  d  i  s  p  ( '  I  e  n  g  t  h  n  o  r  ma  I  i  z  a  t  i  o  n 1  ) 

%  for  icnt  cols  =  l:size(phi,2); 

%  p"Tii(:,icnt  cols)  = 

phi  (:,  icnt  cols),  /nor  m(  phi  (:,  icnt  cols), 'Fro'); 

%  "  end 


end 


END  fNormal  i  ze.  m 


1)))); 
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N.  FOSET  FROM  ASET 


%  %  ********************  fOset  from  Aset.m  *********************** 
% 

%  function  [oset]  =  fOset  from  Aset(ndof,aset); 

%  % 

%%  Usage:  [oset]  =  fOset  from  Aset(ndof,aset); 

%  % 

%  %  T  h  i  s  function  d  e  t  e  r  mi  n  e  s  the  c  o  mp  I  e  me  n  t  a  r  y  subset  "oset" 
%%froma  set  [ 1 :  1 :  n  d  o  f ]  and  the  subset  aset  =  [x  x  x  ...]. 

%  % 

%%  ndof:  Total  number  of  DOF.  Set  is  labeled  "nset". 

%%  aset:  Retained  DOF  (proper  subset  of  [  1 :  1 : n d of ] ) 

%  %  oset:  aset  U  oset  =  n 
%  % 

%  %  Pr  ovi  ded  by  Prof  Gor  di  s 

%  % 

% 

%  n  s  e  t  =  [  1 :  n  d  o  f  ] ; 

% 

%  f  o  r  i  c  n  t  =  1  :  I  e  n  g  t  h  ( a  s  e  t ) ; 

%  indices(icnt)  =find(nset  ==  aset(icnt)); 

%  end 
% 

%  bool  =  ones(si  ze(nset)); 

%  bool(indices)  =  zerosjsize(indices)); 

%  oset  =  nset ( f i  nd( bool  >0) ) ; 

% 

%  %  ********************  fOset  from  Aset.m  *********************** 
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O.  FSPRINGMASS2 


%  7°  fSpri  ngMass2.  m 

%  f  u n c t  i  on  [  k ,  m]  =f  S p r  i  n g Ma s  s  2 ( s  p r  i  n g s ,  mass,  B C ) ; 

%  % 

%%  Usage:  function  [k,  m]=fSpri  ngMass2(spri  ngs,  mass,  BC) 

%  % 

%  %  This  function  script  a  s  s  e  mb  I  e  s  the  stiffness  [  k  ]  and  ma  s  s 

%  %  [  m]  ma  t  r  i  c  e  s  for  an  a  s  s  e  mb  I  a  g  e  of  springs. 

%  % 


%  %  A  linear  chain  of  springs  and  ma s s e s  is  a s s u me d . 

%  %  The  number  of  springs  is  defined  by  the  length  of  the  vector 

'springs' 

%  %  and  their  values  by  the  elements  of  'springs.' 

%  % 

%  %  The  number  of  masses  is  defined  by  the  length  of  the  vector  'mass' 

%  %  and  their  values  by  the  e I  e me n t  s  of  '  ma s  s '  . 

%  %  NOTE:  The  number  of  masses  must  equal  to  the  final  number  of 

act i  ve 

%  %  DOF  ( i .  e. ,  af  t  er  BC' s  appl i  ed) . 

%  % 

%  %  Boundary  conditions  are  specified  by  the  vector  'BC.'  This  vector 

%%  contains  the  DOF  numbers  which  are  to  be  restrained. 

%  % 

%%  For  example,  to  build  the  following  system: 

%  % 

%  %  .01  .02  .  015 

%  %  |  --////■■[  m]  ■■////■■[  m]  ■■////■■[  m] 

%  %  5  6  3.4 

%  % 

%  s  p  r  i  n  g  s  =  [  1  1  ] ; 

%  ma  s  s  =[.01.01]; 

%  BC  =  [ 1] 

%  % 


%  %  BEGI  N  SCRI  PT 

%  %  . 

%  % 

% 

%  if  length(  mass)  ==  (length(springs)+l)  -  length(BC); 

% 

%  k  =  zeros( I  ength( spri  ngs) +1, I  ength(spri  ngs) +1) ; 

%  m  =  z  e  r  o  s  ( I  e  n  g  t  h  (  ma  s  s ) ) ; 

% 

%  %  a  s  s  e  mb  I  e  stiffness  ma  t  r  i  x : 

% 

%  r  ows  =  [  0  1] ; 

%  for  i  s  p  r i  n  g  =  1  :  length(springs); 

% 

%  rows  =  rows  +  1; 

% 

%  addthis  =  [spri  ngs(i  spri  ng)  -spri  ngs( i spri  ng) ; -spri  ngs( i  spri  ng) 

spri  ngs ( i spri  ng) ] ; 

%  k  ( rows,  rows)  =  k(  r  ows ,  r  ows )  +  addthis; 

%  end 

% 

%  i  f  ~i  sempt  y(  BC) ; 

%  keep  =  fOset  from  Aset(  I  ength(  spri  ngs)  +1,  BC) ; 

%  k=k(keep,keep);~ 

%  end 

% 

%  %  a  s  s  e  mb  I  e  ma  s  s  ma  t  r  i  x : 

% 

%  m  =  d  i  a  g  (  ma  s  s ) ; 

% 

%  else 
% 

%  di  sp( 1  Error  in  fSpri  ngmass2.  Check  #  masses,  springs,  and  B  C "  s .  '  ) 

%  return 

% 

%  end 

%  ********************  END  fSpri  n g M a s s 2  m  *********************** 

% 
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P.  KINETICENERGY 


%  This  code  solves  for  the  Kinetic  Energy  function  for  structural  systems 
cl  c 

f  o  r  ma  t  long 

%G  i  v  e  n  data 
%  E  =  1 0  e  6 ; 

%  I  1  =  0.  0  1  6  0  6  7  0  2  1; 

%  n  o  mi  n  a  I  area 
%EI  ement  Lenght 
%  r  h  o 

L  =  total  _l  ength/num_el  ements; 

%  Making  my  eigenvector  matrix  a  42  X  number  of  modes  required 

a  =  zer  os (  4  2,  4  0  ) ; 

a(  3:  42,  :  )  =  p h i  _ ba s e ( :  ,  :  )  ; 

%  For  loops  generated  to  select  the  corresponding  eigenvectors  (magnitudes 
%  and  displace  me  nts)  out  of  the  phi  base  matrix  and  then  multiplied  them 
%  wi  t  h  the  corresponding  hermitian  shape  function  second  derivatives, 
for  i  =  l:num  elements 
for  k  =  1 :  n  d  o  f -  2 
di  spl  ef  t  ( i  ,  k)  =  a( 2*i  ■  1,  k) ; 
rotl  ef  t  ( i  ,  k)  =  a  ( 2* i ,  k); 
di  spri  ght(i  ,  k)  =  a ( 2  *  i  +  1,  k) ; 
r  ot  r  i  g ht ( i  ,  k)  =  a  ( 2*  i  +  2,  k) ; 

%  A  f  t  e  r  integrating  the  original  Hermitians  Shape  functions  squared,  I 
%o  b  t  a  i  n  j  p  r  i  me : 

j  p  r  i  me  ( i  ,  k )  =  l/7*(l/LA2*rotleft(i,k)+2/LA3*displeft(i,k)- 
2/LA3*di  spri  ght(i  ,  k)+l/LA2*rotri  ght(i  ,  k))A2*LA7+l/3*(- 
3  /  L  A2  *  d  i  spl  eftji  ,  k)  +3/  L  A2* d i  spri  g  ht  ( i  ,  k )  -  2  /  L*  r  ot  I  eftji  ,  k )  - 

l/L*rotri  ght  j  i  ,  k) ) * (  1/  L  A2*  r  ot  I  ef  t  ( i  ,  k)  +2/  LA3*di  spl  ef  t  ( i  ,  k)- 

2  /  L  A3* d i  spri  g  ht  ( i  ,  k)+l/LA2*rotri  g  ht  ( i  ,  k) )  *LA6+1/  5*(  2*rot  I  ef  t  ( i  ,  k)*(l/LA2*rotl  ef  t  ( i  , 

k)  +2/  LA3*di  spl  eftji  ,  k)  -  2  /  L A3* d i  spri  g  ht  ( i  ,  k)  +1/  LA2*rotri  g  ht  (  i  ,  k))+(- 

3  /  L A2* di  spl  ef  t  ( i  ,  k)  +3/  L  A2* d i  spri  ght  ( i  ,  k)  -  2  /  L*  r  ot  I  eftji  ,  k)  - 

l/L*rotri  ght(i  ,  k))A2)*LA5+l/4*(2*di  spl  eft(i  ,  k)*(l/LA2*rotl  eft(i  ,  k)+2/LA3*di  spl  eft(i 
,  k)  -  2/  LA3*di  spri  ght  ( i  ,  k )  +1/  L  A2*  r  ot  r  i  ght  j  i  ,  k ) )  +2  *  r  o  1 1  ef  t  ( i  ,  k )  *  ( - 

3  /  L  A2  *  d  i  spl  eftji  ,  k) +3  /  LA2*di  spri  ght ( i ,  k)  -  2  /  L*  r  ot  I  ef  t ( i ,  k)- 

1/  L*  r  ot  r  i  g  ht  ( i  ,  k)))*LA4+l/3*(2*di  spl  eftji  ,  k )  *  ( - 

3  /  L  A2  *  d  i  spl  eftji  ,  k)  +3/  L  A2* d i  spri  g  ht  ( i  ,  k )  -  2  /  L*  r  ot  I  ef  t  ( i  ,  k )  - 

l/L*rotri  ght(i  ,  k))+rotl  eft(i  ,  k)A2)*LA3+di  spl  eft(i  ,  k)*rotl  eft(i  ,  k)*LA2+di  spl  eft(i  ,  k) 
A  2  *  L ; 


f  =  lam  base1; 
kin  I  a  m(  i  ,  k)  =  f  ( :  ,  k) ; 
ki  n_l  amsquared  =  ( ki  n_l  am)  .  A2; 

end 

end 

T  =  ( nomi  nal  _area*rho*ki  n_l  amsquared  /  2) .  *j  pri  me 
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Q.  KINETICENERGYABC 


%%%%K  i  n  e  t  i  c  Energy  calculation  of  the  base  beam,  with  one  ABC  applied,  not 
sitffnes/mass  perturbations. 

cl  c 

f  o  r  ma  t  long 

%G  i  v  e  n  data 
%E  =  1 0  e  6 ; 

%l  1  =  0.  0  1  6  0  6  7  0  2  1; 

L  =  total_l  ength/num_el  ements; 

%Ma king  my  eigenvector  matrix  a  42  X  number  of  modes  required 
aa  =  zeros) 42,  39) ; 

Keep  =  [  3 :  pi  n  n  e  d  -  1, pinned  +  1:42];  %T  his  for  20  elements.  The  42  (ndof)  is  to  be 
changed  if  #  elements  changes 
a  a ( Keep, : )  =  phi  a  OS  E  T ( : , : ) ; 

%F o r  loops  generated  to  select  the  corresponding  eigenvectors  (magnitudes 
%and  displace  me  nts)  out  of  the  phi  base  matrix  and  then  multiplied  them 
%with  the  corresponding  hermitian  shape  function  second  derivatives, 
for  i  =  l:num  elements 
for  k  =  1 :  ndof- 3 
di  spl  ef  t  ( i  ,  k)  =  a  a  ( 2  *  i  -  1,  k) ; 

r  ot  I  ef  t  ( i  ,  k)  =  a  a  ( 2 *  i  ,  k ) ; 

di  spri  ght(i  ,  k)  =  a  a  ( 2  *  i  +  1,  k) ; 

rotri gnt(i  ,  k)  =  a  a ( 2  * i  +  2 ,  k ) ; 


%"  i  a  b  c  p  r  i  me 11  corresponds  to  the  integral  part  of  the  strain  energy  equation 

%a  I  ready  calculated  in  a  polynomial  form,  with  the  corresponding  hermitian  second 

der i  vat i  ve 


%t  i  me  s  its  corresponding  eigenvector 
j  abcpri  me(i  ,  k)  =  l/7*(l/LA2*rotleft( 

2  /  L A3* di  spri  ght(i  ,  k)+l/LA2*rotri  ght  ( 

3  /  LA2*di  spl  eft(i  ,  k)  +3/  LA2*di  spri  ght  ( 
1/  L*  r  ot  r  i  ghtfi  ,  k))*(l/LA2*rotl  e  f  t  ( i 


product  indicated. 
,k)+2/LA3*displeft(i 
,  k)  )  A2*  L  A7  +1/  3*(- 
,  k)  -  2/  L*rot I  ef t ( i ,  k) 
k)+2/LA3*di  spl  eft(i  ,  k) 


k)- 


2/LA3*di  spri  ght(i  ,  k)  +l/LA2*rotri  ght(i  ,  k))*LA6+l/5*(2*rotl  eft(i  ,  k)*(l/LA2*rotl  eft(i  , 
k)  +2/  L  A3* d i  spl  eftfi  ,  k)  -  2/  LA3*di  spri  ght(i  ,  k)+l/LA2*rotri  ght(i  ,  k))+(- 
3  /  LA2*di  spl  eft(i  ,  k)  +3/  LA2*di  spri  ghtfi  ,  k) -  2/ L* rot  I  eftfi  ,  k)  - 

l/L*rotri  ght(i  ,  k))A2)*LA5+l/4*(2*di  spl  eftfi  ,  k)*(l/LA2*rotl  eft(i  ,  k)+2/LA3*di  spl  eftfi 


,  k) -  2  /  L  A3* d i  spri  ght(i  ,  k )  +1/  L  A2*  r  ot  r  i  ght(i  ,  k ) )  +2  *  r  o  1 1  eft(i  ,  k)  *  ( - 

3  /  L  A2  *  d  i  spl  eftfi  ,  k)  +3/  L  A2* d i  spri  ghtfi  ,  k )  -  2  /  L*  r  ot  I  eftfi  ,  k)- 

1/  L*  r  ot  r  i  ght  (  i  ,  k)))*LA4+l/3*(2*di  spl  eftfi  ,  k)*(- 
3  /  L  A2* d i  spl  eftfi  ,  k) +3/ LA2*di  spri  ghtfi  ,  k)  -  2/  L*rot  I  eftfi  ,  k )  - 

l/L*rotri  ghtfi  ,  k) )  +r  ot  I  eftfi  ,  k)  A2) *  L A3+di  spl  eftfi  ,  k)  *  r  ot  I  eftfi  ,  k )  *  L  A2  +d  i  spl  eftfi  ,  k) 


A2  *  L ; 


f  =  I  a maOSET1  ; 
kin  I  a  m(  i  ,  k )  =  f  ( :  ,  k ) ; 
ki  n-l  amsquared  =  (kin  I  am)  .  A2 ; 


end 

end 

%  Assembling  the  strain  energy  equation: 

Tabc  =  (nominal  area*rho*kin  I  amsquared/2) .  *j  abcpri  me 
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MULTIABCS  MODE  ERROR 


%T h e  purpose  of  this  program  is  to  graph  the  error  prediction  wrt  number  of 
%mo des  retained.  This  is  done  for  both  base  system  and  the  identified  ABC. 

%T h i  s  is  purely  a  plotting  code 


norm  T=zeros(  78,  20)  ; 
for  t  =1:  78; 

norm  T(t,  1:20)=  T  sens  t  ot  ( t ,  :  )  /  nor  m(  T  sens  tot(t,  :  ),  i  nf); 
end 

x  =  1:20; 
figure(lOl) 

title('Error  Prediction  for  Base  system1);  hold  on 
s ubpl  ot ( 5, 1,  1) 

%  bar(x, (p(: ,  1)),  .  25,  1  r1  );  hoi  d  on 
b  a  r  (  x ,  ( d  v  calculated  basel(:,l)),.25,'b,);hold  on 
b  a  r ( x,  ( d  v ” c  a  I  cul  at  ed“base2( : , 1) ) , .  25, 1  g1 ) ;  hoi  d  on 
%bar(x,(pabc3(:,l)),725,'g'  ) ;  hoi  d  on 

pi  ot  (  (  BC(  1,  3)  /  2)  ,  0,  1  r  Al  ,  (  BC(  1,  3)  /  2)  ,  0,  1  r  h1  ,  (  BC(  1,  3)  /  2)  ,  0,  1  r*1  )  ;  hoi  d  on 

pi  ot  (  (  BC(  2,  3)  /  2)  ,  0,  1  r  Al  ,  (  BC(  2,  3)/  2)  ,  0,  1  r  h1  ,  (  BC(  2,  3)  /  2)  ,  0,  1  r  * 1  )  ;  hoi  d  on 

%pl  ot ( ( BC( 3,  3)/  2)  ,  0,  1  g*1  )  ;  hoi  d  on 
x  I  i  m(  [  0  21]) 
y  I  i  m(  [  - .  1  .1]) 
grid  on 

s  t  e  m(  mass  Ibis,  dv  mass,  1  g 1  ,  'filled');  hold  on; 

st  em(  mass'l  bl  s,  dv“  mas  s ,  1  k1  ) ,  hold  on; 

stem(EI  Ibis,  dv  Ef/r'/filled'Jlhold  on; 
stem(EI"ibls,  d  v  “  E I  ,  '  k 1  ) ;  hold  on 
t  i  1 1  e(  1  MODE  r  )  " 

%  plot(((BC(:,3)-l)/2),0,lgAI,((BC(:,3)-l)/2),0,'ghl,((BC(:,3)- 

l)/2),  0,  1  g*1  ) 

s  ubpl  ot ( 5, 1,  2) 

%  bar(x, (p(:  ,  2)),  .  25,  1  r1  );  hoi  d  on 
b a r ( x , ( d v  cal  cul  ated_basel(:  ,  2)),  .  25,  1  b1  );  hoi  d  on 
barj  x, ( d  v “ c  a  I  cul  at  ed~base2( : , 2) ) , .  25, 1  g1 ) ;  hoi  d  on 
%bar ( x, ( pabc3( : , 2) ) , 725, 1  g1 ) ;  hoi  d  on 

pi  ot  (  (  BC(  1,  3)  /  2)  ,  0,  1  r  Al  ,  (  BC(  1,  3)/  2)  ,  0,  1  r  h1  ,  (  BC(  1,  3)  /  2)  ,  0,  1  r*'  )  ;  hoi  d  on 

pi  ot  (  (  BC(  2,  3)  /  2)  ,  0,  1  r  Al  ,  (  BC(  2,  3)  /  2)  ,  0,  1  r  h1  ,  (  BC(  2,  3)  /  2)  ,  0,  1  r  * 1  )  ;  hoi  d  on 

%pl  ot ( ( BC( 3,  3)/  2)  ,  0,  1  g*1  )  ;  hoi  d  on 

x  I  i  m(  [  0  21]) 
y  I  i  m(  [  - .  1  .1]) 
grid  on 

s  t  e  m(  mass  Ibis,  dv  ma  s  s ,  1  g 1  ,  'filled');  hold  on; 
stem(mass"lbls,  dv“  mas  s ,  1  k1  ) ,  hold  on; 
stem(EI  Ibis,  dv  Ef/r'/filled'llhold  on; 
st  em(  El  "I  bl  s,  d v ” E I  ,  1  k 1  ) ;  hoi  d  on 
t  i  1 1  e ( '  MODES  1:  2r) 

subpl  ot  (  5,  1,  3) 

%  bar(x, (p(: ,  3)), .  25, 1  r1 );  hoi  d  on 
b  a  r  (  x ,  ( d  v  calculated  b  a  s  e  1  (  :  ,  3) ) ,  .  25,  1  b1  ) ;  hoi  d  on 
b  a  r  (  x ,  (dv'cal  cul  ated~base2(: ,  3)), .  25, 1  g1 );  hoi  d  on 
%b  a  r  (  x,  ( pa  be  3  (  :  , 3  ) ) , 7  2  5, 1  g1 ) ;  hoi  d  on 

pi  ot  (  (  BC(  1,  3)  /  2)  ,  0,  1  r  Al  ,  (  BC(  1,  3)  /  2)  ,  0,  1  r  h1  ,  (  BC(  1,  3)  /  2)  ,  0,  1  r*1  )  ;  hoi  d  on 

pi  ot  (  (  BC(  2,  3)  /  2)  ,  0,  1  r  Al  ,  (  BC(  2,  3)/  2)  ,  0,  1  r  h1  ,  (  BC(  2,  3)  /  2)  ,  0,  1  r  * 1  )  ;  hoi  d  on 

%pl  ot  (  ( BC( 3,  3)/  2)  ,  0,  1  g*1  )  ;  hoi  d  on 

x  I  i  m(  [  0  21]) 
y  I  i  m(  [  - .  1  .1]) 
grid  on 

s  t  e  m(  mass  Ibis,  dv  ma  s  s ,  1  g 1  ,  'filled');  hold  on; 

st  em(  mass'l  bl  s,  d  v“  mas  s ,  1  k1  ) ,  hold  on; 

stem(EI  Ibis,  dv  Ef/rVfilled'lihold  on; 
stem(EI"ibls,  d  v  “  E I  ,  '  k 1  ) ;  hold  on 
t  i  1 1  e( 1  MODES  1:  3r) 

subpl  ot  (  5,  1,4) 

%  barfx, (p( : , 4)), .  25, 1  r1 );  hoi  d  on 
b  a  r  (  x ,  ( d  v  calculated  basel(:,4)),.25,'bl);hold  on 
b  a  r ( x ,  ( d  v “ c  a  I  cul  at  ed“base2( : , 4)), .  25, 1  g1 );  hoi  d  on 
%  x,(pabc3(:,4)),.257'gl);hold  on 

pi  ot  (  (  BC(  1,  3)  /  2)  ,  0,  1  r  Al  ,  (  BC(  1,  3)  /  2)  ,  0,  1  r  h1  ,  (  BC(  1,  3)  /  2)  ,  0,  1  r  * 1  )  ;  hoi  d  on 

pi  ot  (  (  BC(  2,  3)  /  2)  ,  0,  1  r  Al  ,  (  BC(  2,  3)  /  2)  ,  0,  1  r  h1  ,  (  BC(  2,  3)  /  2)  ,  0,  1  r*1  )  ;  hoi  d  on 
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%pl  ot ( ( BC( 3,  3)/  2)  ,  0,  1  g*1  )  ;  hoi  d  on 
x  I  i  m(  [  0  21]) 
y  I  i  m(  [  ■ .  1  .1]) 
grid  on 

s  t  e  m(  mass  Ibis,  dv  mass,  1  g 1  ,  'filled');  hold  on; 
st  em(  mass"l  bl  s,  d  v“  mas  s ,  1  k1  ) ,  hold  on; 
stemj  El  Ibis,  dv  ElVr'.'filled'jjhold  on; 
st  em(  El  "I  bl  s,  d  v  ”  E I  ,  '  k 1  ) ;  hold  on 
t  i  1 1  e ( '  MODES  1: 4r) 

s  ubpl  ot ( 5, 1,  5) 

%  bar(x, (p( : ,  5)),  .  25,  1  r1  );  hoi  d  on 
b  a  r  (  x ,  ( d  v  calculated  basel(:,5)),.25,'bl);hold  on 
barj  x, ( d  v ” c  a  I  cul  at  ed"base2( : ,  5) ) , .  25, 1  g1 ) ;  hoi  d  on 

%b a  r  (  x,  ( pa  be  3  (  :  ,  5  ) ) , 7  2  5, 1  g1 ) ;  hoi  d  on 

pi  ot  (  (  BC(  1,  3)  /  2)  ,  0,  1  r  Al  ,  (  BC(  1,  3)  /  2)  ,  0,  1  r  h1  ,  (  BC(  1,  3)  /  2)  ,  0,  1  r*1  )  ;  hoi  d  on 

pi  ot  (  (  BC(  2,  3)  /  2)  ,  0,  1  r  Al  ,  (  BC(  2,  3)  /  2)  ,  0,  1  r  h1  ,  (  BC(  2,  3)  /  2)  ,  0,  1  r*1  )  ;  hoi  d  on 

%pl  ot ( ( BC( 3,  3)/  2)  ,  0,  1  g*1  )  ;  hoi  d  on 
x  I  i  m(  [  0  21]) 
y  I  i  m(  [  - .  1  .1]) 
grid  on 

s  t  e  m(  mass  Ibis,  dv  ma  s  s ,  1  g 1  ,  'filled');  hold  on; 
st  em(  mass"l  bl  s,  dv"  mas  s ,  1  k1  ) ,  hold  on; 
s  t  e  m(  El  Ibis,  dv  ElVr'.'filled'jjhold  on; 
st  em(  El  "I  bl  s,  d v “ E I  ,  1  k 1  ) ;  hoi  d  on 
t  i  1 1  e( 1  MODES  1:  5r) 


%  %  Base  and  ABC  syst em/ //////////////////// yo  lo 
coment e /////////////////////////////  / 

% 


x  =  1:  20 
f i  gur  e(  103) 

ba r ( x,  ( bbbbf : , 9) ) , .  25, 1  b1 ) ;  hoi  d  on 
%  b a r ( x , ( c c c c ( : , 6 ) ) , .  2 5 ,  1  g 1  ) ;  h o I  d  on 
pi  ot  (  (  BC(  1,  3)  /  2)  ,  0,  1  r  Al  ,  (  BC(  1,  3)/  2)  ,  0, 
pi  ot  (  (  BC(  2,  3)  /  2)  ,  0,  1  r  Al  ,  (  BC(  2,  3)  /  2)  ,  0, 

pi  ot  (  (  BC(  3,  3)  /  2)  ,  0,  1  rAI  ,  (  BC(  3,  3)  /  2)  ,  0, 

pi  ot  (  (  BC(  4,  3)  /  2)  ,  0,  1  r  Al  ,  (  BC(  4,  3)/  2)  ,  0, 

pi  ot  (  (  BC(  5,  3)  /  2)  ,  0,  1  r  Al  ,  (  BC(  5,  3)  /  2)  ,  0, 

grid  on 
x  I  i  m(  [  0  21]) 
ylim([0  .1]) 

s t  e m(  mass  Ibis,  d v  ma s s ,  1  y 1  ,  1  f  i  I  I  e d 1  ) ; 
st  em(  mass"l  bl  s,  d  v“  mas  s ,  1  k1  ) ,  hold  on; 
stemj  El  Ibis,  dv  ElVc'.'filled'jjhold 
st  em(  El  "I  bl  s,  d  v  ”  E I  ,  1  k 1  ) ;  hold  on 
t i 1 1  e( 1  FI VE  ABC  SYSTEMS  MODES  1:  3'  ) 

% 


r  h 1  ,  (  BC(  1,  3)12)  ,  0, 
r  h 1  ,  (  BC(  2,  3)  /  2)  ,  0, 
r  h 1  ,  (  BC(  3 ,  3 ) / 2 ) ,  0 , 
rh1  ,  (  BC(  4 ,  3)  /  2)  ,  0, 
r  h 1  ,  (  BC(  5,  3)  /  2)  ,  0, 


hoi  d  on; 
on; 


r  * 
r  * 
r  * 
r  * 
r  * 


) ;  hoi  d  on 
) ; h o I  d  on 
) ;  h o I d  on 
) ; h o I  d  on 
) ; h o I  d  on 


150 


s. 


NORMRUNTHRU  CRS 


% 


nor  mRUNt  hr  u_c  r  s .  m 


%  This  program  finds  the  NORM  of  the  columns  of  sensitivity  matrix  and  the 
%  NORM  of  the  rows  of  the  inverse  of  the  sensitivity  matrix.  This  was  used 
%  t  o  find  a  correlation  of  the  good  prediction  to  the  ABC  system  used.  It 
%  also  plots  the  information  in  helpful  graphes. 

% 

%  This  program  was  written  for  a  system  of  19  natural  freq.  Set  of  5  modes 
%were  used  in  each  ABC  system,  i.e.,,  modes  1-5,  modes  6-10, or  modes  11-15. 
%  This  accounts  for  the  3  sets  of  modes  per  condition  as  listed  below  in 
%  the  "for"  loop  modeN  =  1:3.  This  program  comp  ares  the  use  of  the  first  5 
%  modes  of  the  base  system  and  one  set  of  five  modes  of  the  ABC  to  that  of 
%  o f  just  10  modes  of  the  ABC.  Notice  that  the  sensitivity  is  a  10x10 
%  square  matrix  indicating  that  only  mass  or  El  changes,  not  both  were 
%  made. 

%  This  program  is  not  part  of  Bui  I  d2Beams  c  r  s .  m  program.  It  is  run 
%  s  e  p  a  r  a  t  e  I  y . 

%  Written  by  Constance  R  S  Fernandez,  Spring  2004 

%  Inputs 

% . 

%  c  o  n  d  b  a  s  e  P I  u  s 
%  i  c  n  t "  o  s  e  t 
%  T  sens  tot 
%  E F  Ibis,  mass  Ibis 
%dv~mass,  d v _ E F ,  dv_cal_ABCten 

%  Out  put  s 

% . 

%  BASE 
%  BASET 

%  abcN,  count  N,  a  cN 
%  modeN,  start  modeN,  start  mode  NT 
%  bb,  t,  tl  NV,  cc,  T,  T I  NV,  vv,  tt 
%  model  a  be  I  NORM 

%normvectT,  normvecTABC,  norm  vecTinvABC 
%  n  o  r  mC 

%  baseN,  baseABCN,  abc  conN,  abc  conTN 
%  baseABCNten,  abc_conNten,  abc_conTNten 

% - Start  program - 


BASE  =  int2str(cond_basePlus(l));  %  cond  no.  of  the  base  line  system 
BASET  =  s  p  r  i  n  t  f  (  1  Base  [1:5]  Cond  =  %s 1  ,  BASE);  %  used  for  plotting 


%  ======|  ni  t  i  a  I  i  z  at  i  on  =====% 


a  be  N  =  0; 
count  N  =0; 


%  =======Ca I  cul  at i  ons  of  NORM  vectors  ======% 

for  count  =  1 :  i  c  n  t  oset  +1  %  number  of  conditions  (base  +  ABC) 
a  _  c  N  =  1; 

for  modeN  =  1:3  %3  sets  of  modes  per  boundry  condition 

start  mo  d  e  N  =  abcN  +  a  _  c  N ;  %  the  beginning  mode  number  of  each  set 

%  indicates  the  use  of  the  first  5  modes  of  base  system  +  5  modes 
%  of  ABC  system 

bb  =  [1:5,  start  modeN:  start  modeN +4]; 

%  labeling  of  modes  for  plotting 
model  abel  NORM  =  i  nt  2s t  r  ( a_ c  N:  a_ c  N+4) ; 

a  _  c  N  =  a  _  c  N  +  5 ;  %a  d  v  a  n  c  e  s  to  the  next  set  of  modes 

%  base  system  plus  5  modes  of  ABC 

t  =  T  sens  t  o  t  (  b  b ,  :  ) ;  %  10x10  matrix 

1 1  NV  =  i  nv[T_sens_t ot  (  bb,  :  ) ) ;  %  10x10  matrix 

%  f  i  r  s  t  10  modes  of  ABC  solo 

st  art  mode  NT  =  abcN  +  1;  %  t  he  beginning  mode  number  of  each  set 
cc  =  [  start  modeNT:  st  a  r  t  modeNT+9] ; 
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T  =  T  sens  tot(cc,:);%  10x10  ma  t  r  i  x 
T I  N  V  =  i  nv[T_sens_t ot  ( cc,  :  ) ) ;  %  10x10  matrix 

%  f  or  loop  for  NORM  of  col  umns  and  r  ows  of  i  nv(  sens  mat  r  i  x) 
for  vv  =  1:10  %  10  rows,  10  columns 
%  base  +  ABC  system 

norm  vecT(  vv,  count  N  +  modeN)  =  n  o  r  m(  t  (  :  ,  v  v ) ) ;  %  columns 
nor  nfvecTi  nv(  vv,  count  N  +  modeN)  =  nor  m(  1 1  NV(  vv,  :  ) ) ;  %  rows 
%  for  10  modes  ABC  solo 

norm  vecTABCf  vv,  count  N  +  modeN)  =  norm(T(:,vv));  %  columns 
nor  m“vecTi  nvABC(  vv,  count  N+modeN)  =  no  r  m(  Tl  N  V  ( v  v ,  :  ) ) ;  %  rows 

end  %  vv  loop 
end  %  ModeN  loop 

abcN  =  abcN+19;  %  advances  to  the  next  ABC  system 
count  N  =  count  N  +  3;  %  counts  up  each  set  of  ABC 
end  %  count  loop 
nor  mC  =4 ;  %  initialize  for  plots 

% . PLOTTI  NG . % 

for  1 1  = 1 : 10  %  figures  (30-40)  plots  6  graphes  per  figure 
f  i  g  u  r  e  ( 1 1  +3  0 ) 
s  u  b  p  I  o  t  ( 3 ,  2 ,  1 ) 
barf  norm  vecT(  :  ,  normC) ) 
title  ('Norm  col  Tsens,  ABC  Modes  [1:5]'); 

subpl  ot  ( 3,  2,  3) 

bar(  norm_vecTi  nv( :  ,  normC) ) 

title  ('Norm  row  TsensINV,  ABC  Modes  [1:5]'); 

subpl  ot  ( 3,  2,  2) 

b  a  r  ( norm  vecTABC(  :  ,  normC) ) 

title  ('Norm  col  Tsens,  ABC  Modes  [1:10]') 

subpl  ot  ( 3,  2,  4) 

b  a  r  ( norm  vecTi  n  v  A  B  C  (  :  ,  normC)) 

title  ('Norm  row  TsensINV,  ABC  Modes  [1:10]') 

subpl  ot  (3,  2,  5) 

%  plotting  error  prediction 

%  using  [1:5]  modes  of  ABC  system  +  [1:5]  modes  of  Base  system; 
baseABCN  =  bar ( d  v  _  c  a  I  _  B  a  s  e  P I  us(  :  ,  nor mC) ,  .  5,  '  r '  ) ;  hoi  d  on 

%  plotting  error  prediction  using  [1:5]  modes  of  Base  system; 
baseN  =  bar ( d v_  c  a  I  _ABC( : , 1) , .  25, '  b' ) ; 


abc  conN  =  int2str(cond  basePI  us(  normC) ) ;  %  cond  no.  for  legend 
abc~conTN  =  spr i  ntf ( '  Basel  1:  5] +ABC[ 1: 5]  Cond  =  %s '  ,  abc_conN); 

grid  on 

I  egend(  [  baseABCN,  baseN] ,  BASET, abc_conTN) ,  hold  on 


%  plotting  actual  error 

if  El  Ibis  ~=[  ]  &  mass  Ibis  ~=[  ] 

s  t  e  m(  mass  Ibis,  d  v  “  ma  s  s ,  '  y '  ,  'filled');  hold  on; 

s  t  e  m(  ma  s  s  ~  I  b  I  s ,  d  v "  ma  s  s ,  '  k '  ) 

Elplot  =  El  I  bl  s  +10;  hoi  d  on 
s  t  e  m(  El  Ibis,  dv  EI,'c','filled');hold  on; 
st  em(  El  ~l  bl  s,  d  v  ”  E I  ,  '  k '  ) 

el  sei  f  mass  Ibis  ~=[  ]  &  El  Ibis  ==[  ] 

s  t  e  m(  ma  s  s  Ibis,  dv  mass7‘  y1  .  1  fi  I  I  ed'  );  hoi  d  on; 

s  t  e  m(  ma  s  s "  I  b  I  s ,  d  v "  ma  s  s ,  '  k '  ) 


else 

s  t  e  m(  El  Ibis,  dv  EI,'c','filled');hold  on; 
stem(  El  ~l  bl  s,  dv~EI  ,  '  k '  ) 

end  %  i  f  E I  _ I  b I  s  ~  =  [  ]  &  ma s s _ I  bl  s  ~=[  ] 

title  ( s  pr  i  nt  f  ( '  Er  r  or ,  Base  [1:5]  +  ABC  [1:5],  pinned  NODE  #  %d '  ,  ( 1 1  +1 ) ) ) 
subpl  ot(3,  2,  6) 

%  plotting  error  prediction  using  [1:10]  modes  of  ABC  system; 
baseABCNten  =  bar ( d  v  _  c  a  I  _  A  B  C t en(  :  ,  nor mC) ) ;  hoi  d  on 
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abc  conNten  =  i  nt2str(cond_ABCten(  normC) ) ; 

abc'conTNten  =  spri  ntf( 1  ABC[  1:  10]  Cond  =  %s 1  ,  abc_conNten) ; 

grid  on 

I  egend( [ bas eABCNt en] , abc_conTNt en) ,  hold  on 

%  plotting  actual  error 

if  El  Ibis  ~=[  ]  &  mass  Ibis  ~=[  ] 

s  t  e  m(  mass  Ibis,  d  v  ~  ma  s  s ,  1  y 1  ,  'filled');  hold  on; 

s  t  e  m(  ma  s  s  ~  I  b  I  s ,  d  v  ~  ma  s  s ,  1  k 1  ) 

El  pi  ot  =  El  lbls+10;hold  on 
s  t  e  m(  El  Ibis,  dv  El/c'.'filled'ljhold  on; 
stem(  El  “I  bl  s,  d v ” E I  ,  1  k 1  ) 

el  sei  f  mass  Ibis  ~=[  ]  &  El  Ibis  ==[  ] 

s  t  e  m(  mass  Ibis,  dv  mass7‘  y1  .  1  fi  I  I  ed1  );  hoi  d  on; 

s  t  e  m(  ma  s  s "  I  b  I  s ,  d  v "  ma  s  s ,  1  k 1  ) 


else 

s  t  e  m(  El  Ibis,  dv  El.'c'.'filled'ljhold  on; 
stem(  El  [I  bl  s,  dv~EI  ,  1  k 1  ) 

end  %  El  _l  bl  s  ~=[  ]  &  ma s s _ I  b I  s  ~=[  ] 


title  ( s  pr  i  nt  f  ( 1  Er  r  or ,  ABC  Modes  [1:10],  pinned  NODE  #  %d 1  ,  ( 1 1  +1 ) ) ) 

normC  =  normC +3;  %  advances  to  the  next  ABC  system 
end  %  1 1  =  1:10  for  plotting  graphes 

o/0  ********************  END  nor  mRUNt  hr  u  crs  m  *********************** 
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T. 


PLOTBEAMMODES  CRS 


Plot  BeamModes_c  r  s .  m 
Cal  cul  at  es  nat  ur al  f  r equenci  es 
Provided  by  Prof  Gor  di  s 
I  nput  s  needed: 


k  a ,  ma , 

mx,  kx 

P  r  o  g  r  a  ms 

needed : 

f  Modes 

Out  put  s : 

%  I  a  ma ,  phi  a,  lamx,  phix  (without  rigid  body  modes) 
%  n  u  m  r  b  m 

%  p  h  i  a  plot,  phix  plot 


di  sp( 1  '  ) ; 

d  i  s  p  ( '  Calculating  modes  for  each  beam.  ..plot  frequency  comparison1) 


%  Get  modes  of  each  beam: 

[  I  ama,  phi  a]  =f  Modes)  ka,  ma) ; 

[  I  amx,  phi  xj  =f  Modes ( kx,  mx) ; 

%used  to  plot  the  mode  shapes  org  BC  before  ABC 
phia  plot  =  phia; 
phi  x  “  p  I  o  t  =  phix; 

%  Set  any  rigid  body  mode  freqs  to  zero: 
num_rbm  =  I  ength(fi  nd(l  ama  <  1)); 

spri  ntf('  Number  of  Rigid  Body  Modes  Found:  %2 i  1  ,  num_rbm) 

disp(  1  Removing  rigid  body  mode  frequencies  from  vectors...1) 

I  a  ma  =  I  a  ma  ( f  i  n  d  ( I  a  ma  >  1 ) ) ; 

I  a  mx  =  I  a  mx  ( f  i  n  d  ( I  a  mx  >  1 ) ) ; 

o/0  ********************  end  PI  ot  BeamModes  crs.m  *********************** 
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U.  PLOTTINGBARS  CRS 


%  ********************  p| otti  ngBARScrs.  m 
%To  be  used  with  Build2Beams.m  and 


%  This  program  plots  9  graphes  per  figure.  The  first  columns  of  3  graphes 
%  are  the  mode  shapes  of  the  ABC  system  used  in  error  prediction.  The  next 
%  c  o  I  u  mn  of  3  graphes  are  the  error  prediction  using  only  5  modes  of  ABC 
%  system.  The  last  column  of  3  graphes  are  the  error  predictions  using  the 
%  first  5  modes  of  base  system  plus  5  modes  of  the  ABC  system.  The  row 
%  represent  modes  1-5,  middle  row:  modes  6-10,  last  row:  modes  11-15.  Each 
%  o f  the  error  prediction  graphes  also  have  the  base  only  prediction  and 
%  t  h  e  actual  error  plotted  for  easy  reference. 

% 

%  Written  by  Constance  R  S  Fernandez,  Spring  2004 

%  Inputs 

% . 

%  c  o  n  d  b  a  s  e  P I  u  s 
%  FOM  A B C 5 p e r ,  FOM  PLUSper 
%  i  c  n  t  o  s  e  t 
%  modes h ape 
%  rel  freqERROR 
%  ypos 

%  E I  Ibis,  mass  Ibis 
%  d  v  ~  ma  s  s ,  d  v  _  E  F 

%  Out  put  s 

% . 

%  BASE,  BASET 

%  FOMBASE,  FOMABC,  FOMPLUS 
%  i  n  t  e  r  v  e  I  p 

%  E R ,  barp,  shape,  error,  a  cp,  modep,  ap, 

%  model  a b e I  p ,  FOMABCI  abel  p,  "FOMPLUSI  abel  p 
%  a  be  con,  a  be  conT 
%  A B C 7  base 
%  plus  con,  plus  conT 
%  base?  plus,  El  “pi  o  t 


BASE  =  int2str(cond  ba s e P I  u s ( 1 ) ) ; 

FOMBASE  =  i  n 1 2  s  t  r  (  FOM  ABC5per(  1)  )  ; 

BASET  =  s  pr  i  nt  f (  1  Bas  e"Cond  =  %s,  FOM  =  %s 1  ,  BASE,  FOMBASE); 


i  nt  er  vel  p  =  3; 
modes h ape  =  1; 

ER  =  1; 

for  barp  =  l:icnt_oset 


f  i  gur e( bar p+10)  %  figures  11-20 
f  o  r  ma  t  bank 

%s  h  a  p  e  =  [  modes  ha  pe:  modes  ha  pe +2  0] ;  %Thi  s  number  is  equal  to  number  of  elements, 
shape  =  [  modes  ha  pe:  modes  ha  pe +20] ;  %Thi  s  number  is  equal  to  number  of  elements, 
error  =  r  o  u  n  d ( r  e I  f  reqERROR( ER:  ER  +  15) *100) /  100; 
a  cp  =  1; 

for  modep  =  1:3  %3  sets  of  modes  per  boundry  condition 


ap  =  [a  cp:  a  cp+4] ;  %modes 
%  REL  errorl  =  int2str(error(a  cp)); 

%  Errorl  abel  p  =  sprintff'Rel  error  =  %s 1  , 

%  REL  err  or  2  =  i  nt 2st r ( er r or ( a  c  p  +1 ) ) ; 

%  Errorl  abel  p  =  s pr  i  nt f ( 1  Rel  error  =  %s 1  , 

%  REL  e  r  r  o  r  3  =  int2str(error(a  c  p  +2 ) ) ; 

%  Errorl  abel  p  =  s pr  i  nt f ( 1  Rel  error  =  %s 1 , 

%  REL  err  or  4  =  i  nt 2st r ( er r or ( a  c  p  +3 ) ) ; 

%  Errorl  abel  p  =  s  p  r  i  n  t  f ( 1  R  e I  error  =  %s 1  , 

%  REL  errors  =  i  nt 2st r ( er r or ( a  cp+4)); 

%  Errorl  abel  p  =  s pr  i  nt f ( 1  Rel  error  =  %s 1  , 


REL_er  r  or  1) ; 
REL_er  r  or  2) 
RE  L_  e  r  r  o  r  3 ) 
REL_er  r  or  4) 
REL_er  r  or  5) 


modelabelp  =  i  n  1 2  s  t  r  ( a  c  p :  a  cp+4); 

FOMABC  =  int2str(FOM  ABC5per(  i  ntervel  p+modep) )  ; 
FOMABCI  abel  p  =  spri  ntf('  System  FOM  =  %s 1  ,  FOMABC); 


FOMPLUS  =  int2str(FOM  PLUSper ( i  nt er vel  p+modep) ) ; 
FOMPLUSI  abel  p  =  spri  ntf('  System  FOM  =  %s 1  ,  FOMPLUS); 
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% 


■% 


=mode  shape  or  beam  X  and  beam  A  with  A B C  = 


f i gur  e( bar  p+50) 
subpl  o t  ( 3 ,  1,  modep) 
pi  ot ( ypos ,  di  s  pA  t  ot ( s  hape, a 
di  s  pA  t  ot ( s  hape,  a  c  p+1) , 
di  s pA~t  ot ( shape, a-cp+2) , 
di  s  pA"t  ot ( s  hape,  a“c  p+3) , 
di  s  pA"t  ot ( s  hape,  a"c  p  +4 ) , 
di  s  pX~t  ot ( s  hape, a"c  p) , 
di  spX“t  ot ( shape,  a"c p+1) , 
di  spX“t  ot ( shape,  a~c  p+2) , 
di  s  pX~t  ot  j  s  hape, a~c  p+3) , 
di  spX“t  ot ( shape,  a"c p+4) , 


k-  o'  ,  y pos , 


g- 5 ; , 

ypos,  ... 

b-  d  , 

ypos,  ... 

r-  x'  , 

ypos,  ... 

m-  * '  , 

ypos,  ... 

r  -  -  o' 

ypos,  .. 

b-  -  s' 

.  ypos,  .. 

m-  -  d ' 

.  ypos,  .. 

c  -  -  x ' 

,  ypos,  .. 

k-  -  * ' 

) ,  gr i d  on 

legend(sprintf('Rel  Freq  Error  =  %d',  e  r  r  o  r  ( a  cp)),... 
sprintff'Rel  Freq  Error  =  %d 1  ,  error)  a  c  p+1 )),... 

sprintfj'Rel  Freq  Error  =  %d'  ,  er r or j a~c p+2) ) , . . . 

sprintfj'Rel  Freq  Error  =  %d'  ,  er r or j a~c p+3) ) , . . . 

sprintfj'Rel  Freq  Error  =  %d 1  ,  error)  a“cp+4)  ) ) ; 

%  I  egend(  s  pr  i  nt  f  ( 1  Bm  X,  Md  %d '  ,  a  cp1),  s  pr  i  nt  f  ( 1  Bm  X,  Md  %d 1  ,  a  c  p  +1 ) ,  .  .  . 

%  s  pr  i  nt  f  ( 1  Bm  X,  Md  %d'  ,  a  cp+2),  s  pr  i  nt  f  (  1  Bm  X,  Md  %d '  ,  a  c  p  +3 ) ,  .  .  . 

%  s pr  i  nt f j  1  Bm  X,  Md  %d '  ,  a"cp+4),  spri  ntf( 1  Base,  Md  %d 1  ,  a  cp),  ... 

%  s pr i  nt f j 1  Bas e,  Md  %d 1 ,  a"cp+l),  s pr i  nt f j 1  Bas e,  Md  %d 1 ,  a" cp+2),  ... 

%  s  pr  i  nt  f  j  1  Bas  e,  Md  %d'  ,  a “ c  p  +3 ) ,  s  pr  i  nt  f  j  1  Bas  e,  Md  %d '  ,  a “ c  p  +4 ) ) 


t  i  1 1  e  ( s  p  r  i  n  t  f  ( 1  Mo  d  e  s  [  %s  ] '  ,  model  a  b  e  I  p ) ) 
axi  s  t i  ght 

% . % 

%  =====  bar  graphes  of  error  solution  using  only  ABC=============% 

%  ===============================================================% 

f  i  gur  e( bar  p  +  10)  %  f i  gur  es  11-20 


subpl  ot(3,  2,  2*modep  -  1) 

abc  con  =  int2str(cond  ABC(  i  ntervel  p+modep) ) ; 

abc"conT  =  s pr i  nt f ( 1  ABC  Cond  =  %s,  FOM  =  %s 1  ,  a b c _ c o n ,  FOMABC); 

ABC  =  bar) d  v  cal  ABC):  ,  i  ntervel  p+modep),  .  5,  1  r1  );  hold  on 
%ABC  =  b  a  r  ( d  v  c  a  F  ABC)  :  ,  1 :  5 ) ,  .  5 ,  1  r 1  ) ;  hold  on  %trying  to  isolate  the  first 
f  i  v  e  mo  d  e  s 

base  =  b  a  r  ( d  v  cal  ABC)  :  ,  1) ,  .  25,  1  b1  ) ;  hoi  d  off  %o  n  %  base  first  5  modes 
%F OMa  be  =  bar[l,0l;  hold  off 
grid  on 

%l  egend)  [  ABC,  base,  F  OMa  be  ]  ,  pi  us  conT,  BASET,  FOMABCI  abel  p)  ,  hold  on 
%  g r i  d  on 

I  egend) [ ABC,  base] , abc_conT,  BASET) ,  hold  on 
ti  tl  e(spri  ntf) 1  ABC  only,  [  %s]',  model  abel  p) ) ; 


if  El  Ibis  ~=0  &  mass  Ibis  ~=0 
%“  p  I  o  t  s  actual  error 

stem)  ma  s  s  Ibis,  dv  mass,  1  y 1  ,  'filled');  hold  on; 
stemjmass'lbls,  dv“mass,'k'),  hold  on; 

Elplot  =  El  Ibis +10;  hold  on  %  last  half  of  plot 
stem)  El  plot  7  dv  EI,'c','filled');hold  on; 
stem)  Elplot,  dv"EI,'k');  hold  on 

%  plots  the  green  triangle  which  indicates  pinned  node 
%pl  ot(barp+9,  0,  '  gAI  ,  barp+9,  0,  '  gh'  ,  barp+9,  0,  '  g*'  ,  .  .  . 

%  bar  p +11,  0,  '  gAI  ,  bar  p  +  1 1 ,  0,  1  gh'  ,  bar  p +  11,  0,  '  g*'  ) 

plot(((BC(:,3)-l)/2),0,'gA',((BC(:,3)-l)/2),0,'gh',((BC(:,3)- 
1) / 2)  ,  0,  1  g*' , .  .  . 

((BC(:,3)-l)/2),0,'gA',((BC(:,3)-l)/2),0,'gh',((BC(:,3)-l)/2),0,'g*' 


el  sei  f  mass  Ibis  ~=0  %&EI  Ibis  =0 
%  plots  actual  error 

stem)  ma  s  s  Ibis,  dv  mass,  '  y 1  ,  'filled');  hold  on 
stemjmass"lbls,  d  v  ~  ma  s  s ,  '  k 1  ) ;  hold  on 
%  plots  the  green  triangle  which  indicates  pinned  node 
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1) / 2) ,  0,  1  g* 


%pl  ot  (  bar  p  +9 ,  0,  1  g Al  ,  bar  p  +9 ,  0,  1  g  h1  ,  bar  p  +9 ,  0,  1  g* 1  ) 
p  ot(((BC(:,3)-l)/2),0f,g^,/((BC(:f3)-l)/2)f0f,gh,,((BC(:,3)- 


el  s e 

%  p I  o  t  s  actual  error 

s  t  e  m(  El  Ibis,  dv  El,  1  c 1  ,  'filled'  ) ;  hoi  d  on; 
s  t  e m(  E I  ~  I  bl  s ,  d  v  “  E I  ,  1  k 1  ) ;  hold  on 

%  plots- the  green  triangle  which  indicates  pinned  node 
%pl  ot  (  bar  p  +9 ,  0,  1  gAI  ,  bar  p  +9 ,  0,  1  gh1  ,  bar  p  +9 ,  0,  1  g*1  ) 
plot(((BC(:,3)-l)/2),0,lgAI,((BC(:,3)-l)/2),0,lghl,((BC(:,3)- 
l)/2),0,'g*1) 

end  %  El  Ibis  ... 


%  ===============================================================% 

%  =========b a r  graphes  of  error  solution  using  ABC  +  ba s e ========% 

% . % 

subpl  o  t  ( 3 ,  2,  2  *  mo  d  e  p ) 


plus  con  =  int2str(cond  basePI  us(i  ntervel  p+modep));  %  for  legend 
plus"conT  =  spri  ntf('  Base+ABC  Cond  =  %s  FOM  =  %s',  plus  con,  FOMPLUS);%  for 

legend 

plus  =  b  a  r  ( d  v  cal  Bas  ePI  us  (:,  i  nt  er  vel  p+modep) ,.  5,  1  r 1  ) ;  hold  on 
base  =  bar  ( dv"cal"ABC(  :  ,  1) ,  .  25,  1  b1  ) ;  hold  on  %  base  first  5  modes 
%F  OMp  I  us  =  b  a  r  (  1 ,  D )  ;  hold  off 
grid  °  n 

I  egend( [ pi  us,  base] ,  pi  us  conT,  BASET) ,  hold  on 
%  I  egend(  [  pi  us,  base,  F  OMp  F us] ,  pi  us_conT,  BASET,  FOMPLUSI  abel  p) ,  hold  on 

if  El  Ibis  ~=0  &  mass  Ibis  ~=0 
%“  p  I  o  t  s  actual  error 

stem(mass  Ibis,  dv  ma  s  s ,  1  y 1  ,  'filled');  hold  on; 
st  em(  mass"l  bl  s,  d v"  ma s  s ,  1  k 1  ) ;  hold  on 
Elplot  =  El  Ibis +10;  hold  on  %  last  half  of  plot 
stem(EI  pi  ot7  dv  El.'c'/filled'ljhold  on; 
st  em(  Elplot,  dv"EI,'kl);hold  on 

%  plots  the  green  triangle  which  indicates  pinned  node 
%pl  o t  ( bar  p  + 1 ,  0,  1  gAI  ,  bar  p  + 1 ,  0,  1  gh1  ,  bar  p  + 1 ,  0,  1  g*1  ,  .  .  . 

%  barp+11,  0,  1  gAI  ,  barp+11,  0,  1  gh1  ,  barp  +  11,  0,  1  g*1  ) 

plot(((BC(:,3)-l)/2),0,lgAI,((BC(:,3)-l)/2),0,lgh',((BC(:,3)- 
1) / 2) ,  0,  1  g*1 , .  .  . 

((BC(:  ,  3)-l)/2),  0,  1  gAI  ,  ((BC(:  ,  3)-l)/2),  0,  1  gh1  ,  ((BC(:  ,  3)-l)/2),  0,  1  g*1 


el  sei  f  mass  Ibis  ~=0  %&EI  Ibis  =0 
%  p I  o  t  s  actual  error 

stem(mass  Ibis,  dv  mass.'y'/filled'lihold  on; 
st  em(  mass"l  bl  s,  dv"mass,'k');hold  on 
%  plots  the  green  triangle  which  indicates  pinned  node 
%pl  ot  (  bar  p  +9 ,  0,  1  g Al  ,  bar  p  +9 ,  0,  1  gh1  ,  bar  p  +9 ,  0,  1  g* 1  ) 
plot(((BC(:,3)-l)/2),0,lgAI,((BC(:,3)-l)/2),0,'ghl,((BC(:,3)- 

1) / 2) ,  0,  1  g*1  ) 
else 

%  p I  o  t  s  actual  error 

stem(EI  Ibis,  dv  El.'c'/filled'lihold  on; 
st  em(  El  "I  bl  s,  d  v  “  E I  ,  1  k 1  ) ,  hold  on 

%  plots" the  green  triangle  which  indicates  pinned  node 
%  pi  ot(  barp+9,  0,  1  gAI  ,  barp+9,  0,  1  gh1  ,  barp+9,  0,  1  g*1  )%changed  bar  p  +1 
pi  o  t  (  (  (  B  C  (  :  ,  3 )  -  1 )  /  2 )  ,  0,  1  gAI  ,  (  (  B  C  (  :  ,  3 )  -  1 )  /  2 )  ,  0,  1  gh1  ,  (  (  B  C  (  :  ,  3 )  - 

1) / 2) ,  0,  1  g*1  ) 

end  %  i  f  E I _ I  b I  s  ... 

title(sprintf('Base  [1:5]  +  ABC  [  %s  ] 1  ,  model  abel  p) ) ; 

%  t  i  1 1  e(  spr  i  ntf  (  1  Base  +  ABC,  pinned  at  NODE#  %d 1  ,  b  a  r  p  +1)) 

a  cp=  a  cp  +5; 
end  %  mo d e p  loop 

mo  deshape  =  mo  deshape  +  11;  %  advances 

intervelp  =  intervelp  +3;  %  advances  to  the  next  ABC  system 
ER  =  ER+19; 
end  %  b  a  r  p  loop 
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%  f  i  gur  e( bar  p  +  11) 

%pl  ot ( 1:  20, ABC1,  bar) 


%  %ABC  Mode  1  Error 
%  subpl  ot(5,  1,  1) 

%  ABC1  =d v  cal  ABC1; 

%  s  t  e  m(  El  F  b  I  s ,  d  v  E I  ,  1  k '  ) ;  hold  on 

%  bar(l:  20,  ABC1,  .  25,  1  r 1  )  ; 

%  plot(((BC(:,3)-l)/2),0,'gAI,((BC(:,3)-l)/2),0,lghl,((BC(:,3)- 

1 )  / 2 ) ,  0,  '  g* 1 ) ; 

%  t i 1 1  e( 1  ABC  Mode  l1  ) 

%  %ABC  Mode  2  Error 
%  subpl  ot  (  5,  1,  2) 

% 

%  ABC2  =  dv  cal  ABC2 ; 

%  s  t  e  m(  El  Ibis,  d  v  E I  ,  1  k 1  ) ;  hold  on 

%  bar(l:  20, ABC2 , .  25, 1  r 1  )  ; 

%  plot(((BC(:,3)-l)/2),0,'gAI,((BC(:,3)-l)/2),0,lgh',((BC(:,3)- 

1 )  /  2 ) ,  0, 1  g* 1  )  ; 

%  title!1  ABC  Mode  V  ) 

% 

%  %ABC  Mode  3  Error 
%  subpl  ot  ( 5,  1,  3) 

% 

%  ABC3  =  dv  cal  ABC3 ; 

%  stem(  El  Ibis,  d  v  E I  ,  '  k 1  ) ;  hold  on 

%  bar(l:  20, ABC3 , .  25, 1  r 1  )  ; 

%  pi  ot(((BC(:  ,  3)-l)/2),  0,  1  gAI  ,  ((BC(:  ,  3)-l)/2),  0,  1  gh'  ,  ((BC(:  ,  3)- 

1) / 2) ,  0,  1  g*1 )  ; 

%  title!1  ABC  Mode  31 ) 

% 

%  %ABC  Mode  4  Error 
%  subpl  ot  ( 5,  1,  4) 

%  ABC4  =d v  cal  ABC4; 

%  stem!  El  Fbls,  d  v  E I  ,  1  k 1  ) ;  hold  on 

%  bar(l:  20, ABC4 , .  25, 1  r 1  )  ; 

%  plot(((BC(:,3)-l)/2),0,'gA,,((BC(:,3)-l)/2),0,lghl,((BC(:,3)- 

1) / 2) ,  0,  1  g*1 ); 

%  title!1  ABC  Mode  41 ) 

% 

%  %ABC  Mode  5  Error 
%  subpl  ot  ( 5,  1,  5) 

%  ABC5  =d v  cal  ABC5 ; 

%  stem!  El  Fbls,  d  v  E I  ,  '  k 1  ) ;  hold  on 

%  bar(l:  20, ABC5 , .  25, 1  r 1  )  ; 

%  plot(((BC(:,3)-l)/2),0,'gA,,((BC(:,3)-l)/2),0,lghl,((BC(:,3)- 

1) / 2) ,  0,  1  g*1 ); 

%  title!1  ABC  Mode  51 ) 


END  pi  ot  t  i  n g B A R S _ c  r  s .  m 
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V. 


STRAINENERGY 


%%%S train  Energy  calculation  of  the  base  beam,  that  is  neither  ABC's  nor 
%%%s itffnes/mass  perturbations  applied. 

cl  c 

f  o  r  ma  t  long 

%G  i  v  e  n  data 
%  E  =  1 0  e  6 ; 

%  I  1  =  0.  0  1  6  0  6  7  0  2  1; 

%EI  ement  Lenght 

L  =  total_l  ength/num_el  ements; 

%Ma king  my  eigenvector  matrix  a  42  X  number  of  modes  required 

a  =  z e r  o s (  4  2,  4  0  ) ; 

a( 3: 42, : )  =  p h i _ ba s e ( : , :  ) ; 

%F o r  loops  generated  to  select  the  corresponding  eigenvectors  (magnitudes 
%and  displace  me  nts)  out  of  the  phi  base  matrix  and  then  multiplied  them 
%with  the  corresponding  hermitian  shape  function  second  derivatives, 
for  i  =  l:num  elements 
for  k  =  1 :  n  d  o  f -  2 
di  spl  ef  t  ( i  ,  k)  =  a( 2*i  -  1,  k) ; 
rotl  eft(i  ,  k)  =  a  ( 2  * i ,  k ) ; 
di  spri  ght(i  ,  k)  =  a ( 2 * i  +  1,  k) ; 
rotri  gnt(i  ,  k)  =  a ( 2  *  i  +  2 ,  k ) ; 

%"  i  p  r  i  me "  corresponds  to  the  integral  part  of  the  strain  energy  equation 

%a  I  ready  calculated  in  a  polynomial  form,  with  the  corresponding  hermitian  second 

der i  vat i  ve 

%t  i  me  s  its  corresponding  eigenvector  product  indicated, 
i  p  r  i  me  ( i  ,  k )  =  l/3*(6/LA2*dispieft(i,k)+2/L*rotleft(i,k)- 

6  /  L A2* di  spri  ght(i  ,  k)  +4/L*rotri  ght  ( i  ,  k ) )  A3/(12/LA3*di  spl  eft(  i  ,  k)  +6/LA2*rotl  eft(i  ,  k)  - 
12/  LA3*di  spri  ght(i  ,  k)+6/LA2*rotri  ght(i  ,  k))-l/3*(-6/LA2*di  spl  eft(i  ,  k)  - 
4  /  L  *  r  o  1 1  ef  t  ( i  ,  k) +6/ LA2*di  spri  ght  ( i ,  kj- 

2/L*rotright(i,k))A3/(12/LA3*displeft(i,k)+6/LA2*rotleft(i,k)- 
12/  LA3*di  spri  ght  ( i  ,  k)  +6/  L  A2*  r  ot  r  i  ght  ( i  ,  k ) ) ; 

s  t  r  a  i  n ( i  ,  k )  =  i  p  r  i  me ( i  ,  k ) ; 

end 

end 

%  Assembling  the  strain  energy  equation: 

UN  =  (nominal  El/2)*strain 
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W.  STRAINENERGYABC 


%%%S train  Energy  calculation  of  the  base  beam,  with  one  ABC  applied,  not 
sitffnes/mass  perturbations. 

cl  c 

f  o  r  ma  t  long 

%G  i  v  e  n  data 
%E  =  1 0  e  6 ; 

%l  1  =  0.  0  1  6  0  6  7  0  2  1; 

L  =  total_l  ength/num_el  ements; 

%Ma king  my  eigenvector  matrix  a  42  X  number  of  modes  required 
aa  =  zeros) 42,  39) ; 

Keep  =  [  3 :  pi  n  n  e  d  -  1, pinned  +  1:42];  %T  his  for  20  elements.  The  42  (ndof)  is  to  be 
changed  if  #  elements  changes 
a  a ( Keep, : )  =  phi  a  OS  E  T ( : , : ) ; 

%F o r  loops  generated  to  select  the  corresponding  eigenvectors  (magnitudes 
%and  displace  me  nts)  out  of  the  phi  base  matrix  and  then  multiplied  them 
%with  the  corresponding  hermitian  shape  function  second  derivatives, 
for  i  =  l:num  elements 
for  k  =  1 :  ndof- 3 
di  spl  ef  t  ( i  ,  k)  =  a  a  ( 2  *  i  -  1,  k) ; 

r  ot  I  ef  t  ( i  ,  k)  =  a  a  ( 2 *  i  ,  k ) ; 

di  spri  ght(i  ,  k)  =  a  a  ( 2  *  i  +  1,  k) ; 

rotri gnt(i  ,  k)  =  a  a ( 2  * i  +  2 ,  k ) ; 

%"  i  a  b  c  p  r  i  me 11  corresponds  to  the  integral  part  of  the  strain  energy  equation 

%a  I  ready  calculated  in  a  polynomial  form,  with  the  corresponding  hermitian  second 

der i  vat i  ve 

%t  i  me  s  its  corresponding  eigenvector  product  indicated, 
i  abcpri  me(i  ,  k)  =  l/3*(6/LA2*displeft(i,k)+2/L*rotleft(i,k)- 

6  /  L  A2* d i  spri  ght(i  ,  k)  +4/  L*rotri  ght(i  ,  k ) )  A3/(12/LA3*di  spl  eft(i  ,  k)  +6/  LA2*rot  I  eft(i  ,  k)  - 
12/  LA3*di  spri  ght  ( i  ,  k)  +6/  LA2*rot  r  i  ght(i  ,  k))-l/3*(-6/LA2*di  spl  e  f  t  ( i  ,  k)  - 
4/  L*  r  ot  I  ef  t  ( i  ,  k)  +6  /  LA2*di  spri  g  ht  ( i  ,  k)  - 

2/L*rotri  ght(i  ,  k))A3/(12/LA3*di  spl  eft(i  ,  k)  +6/  LA2*rot  I  eft(i  ,  k)  - 
12/  LA3*di  spri  ght  ( i  ,  k)  +6/  L  A2*  r  ot  r  i  ght  ( i  ,  k)j; 

s  t  r  a i  n ( i ,  k )  =  i  abcpri  me(i  ,  k); 

end 

end 

%  Assembling  the  strain  energy  equation: 

UNABC  =  (nominal  El/2)*strain; 
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